The Gregoriev Ice Cap evolution according to the 2-D ice flowline model for various climatic scenarios in the future

The Gregoriev Ice Cap evolution according to the 2-D ice flowline model for various climatic scenarios in the future Y. V. Konovalov and O. V. Nagornov Moscow Engineering Physics Institute (State University), Kashirskoe shosse, 31, 115409 Moscow, Russia Received: 27 October 2008 – Accepted: 18 November 2008 – Published: 22 January 2009 Correspondence to: Y. V. Konovalov (yu-v-k@yandex.ru) Published by Copernicus Publications on behalf of the European Geosciences Union.


Introduction
Gregoriev Ice Cap is one of the typical plane-top glaciers in the Central Asia (Figs. 1     and 2).Plane-top glaciers have average ice thickness in the range of 100-150 m.Maximum length and width are equal to 6 and 4 km, respectively.They are located at the South slope of the Terskey Ala-Tau (Fig. 1).The glaciers contain paleoenvironmental records in the Central Asia.Moreover they are sources of fresh water in the region.
Resent studies of the climate variability in the region and the satellite observations revealed that the plane tops glaciers retreated last few decades.So a range of ques- tions arise ahead of glaciological community.Does the retreat respond to a long-term or a shot-term climate changes?Can the plane tops ice caps disappear and what is the time of disappearing?Do the modern environmental circumstances lead to disappearing of the plane tops glaciers or the climate should be warmer?
The Ice Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM, Pattyn et al., 2008) have revealed that modern higher-order (Hindmarsh, 2004) and full Stokes ice sheet models can be successfully applied in the cases of high basal topographic variability and slipperiness when the longitudinal stresses should be taken in to account and the SIA (Hutter, 1983) is not valid.Hence, plane tops glaciers modeling based on the higher-order and full Stokes ice sheet models can assign the glacier retreat/advance correctly and give a forecasts of the glaciers evolution.ISMIP-HOM results have shown that in general all participating models are in concordance.Full Stokes models closely agree with each other.Their spread of results is very low (<1%) and they give a result in concordance with analytical solutions based on perturbation theory for a constant viscosity.In comparison to the full Stokes models the results of the higher-order approximations show larger spread (at the order of 2.5%).The higher spread is found for high aspect rations where all stress components become important.Furthermore, the full Stokes models results anti-correlate with ones of the higher-order models for the highest aspect ratio (0.2, Pattyn et al., 2008).
The evolution of the Gregoriev Ice Cap at this paper was obtained by 2-D ice-sheet modeling and by using the ice surface mass balance measurements which were carried out in 1987and 1988(Fig. 3, Mikhalenko, 1989).The developed 2-D ice sheet model is the full Stokes dynamic flowline model which describes ice flow along a flowline (Pattyn, 2000).It's supposed that the flowline is fixed in space (Fig. 2).The model takes into account divergence/convergence of the ice flow (Pattyn, 2000).
The secondary objective of this study was search of correlations between the glacier length changes and the atmospheric temperature histories similar to equation early obtained by Oerlemans based on in-situ measurements (Oerlemans, 1994)  2 Field equations

Diagnostic equations
Mathematical model for obtaining ice flow velocity in steady-state glacier along the flowline includes the continuity equation for incompressible substance, the mechanical equilibrium equation in terms of stress deviator components and the Glen law (Pattyn, 2000): (1) The description of the constants and variables in the equations of the paper are given in Table 1.The second equation of the system Eq.( 1) stems from the mechanical equilibrium equations for the solid domain under the gravitation ( Finally, taking into account the relations between deviatoric stresses and strain rates, the diagnostic equations can be written in the form of the two integro-differential flow velocity equations:

Conclusions References
Tables Figures

Back Close
Full (2) The flow-law rate factor temperature dependence was taken from Paterson (1994).

Prognostic equation
The prediction of ice thickness changes along the flowline as a response to the time dependent environmental conditions is given by the mass balance equation (prognostic equation, MacAyeal et al., 1996): where ū is the depth averaged horizontal velocity.The last term in Eq. ( 3) provides the stability of the prognostic equation numerical solution.The instability can arises at high spatial resolution and steep velocity decreasing, especially at the ice front vicinity, where ice thickness abruptly vanish.The effect of the last term leads to ice mass disappearance and depends on the artificial viscosity (ν).So it's negative from the point of mass conservation law.And it requires to choose the artificial viscosity value as small as it possible to achieve the numerical stability.
We considered non-zero artificial viscosity values closely to the glacier front only, exactly, at the distance about 200 m upstream from the front point.The harmonic method of the numerical stability analysis leads to the expression for the non-zero viscosity values where ∆x is the step of the spatial grid, u i is the horizontal velocity at the i th grid node and the staggered grid convention is used (MacAyeal et al., 1996).
In fact, the artificial viscosity depends on the velocity gradients at given spatial resolution and increases with the velocity gradients increasing.The ν 0 was assessed in the steady state experiments described below.

Heat transfer equation
The heat transfer in a glacier is a result of the thermal diffusion, horizontal and vertical advection, and deformational heating.The ice temperature changes along a flowline are described by the heat transfer equation (Pattyn, 2000): 3 Boundary conditions 3.1 Boundary conditions at ice surface It's supposed that the boundary conditions for the diagnostic equations Eq. ( 2) also should be written in the form of two flow velocity equations.One of them is the mass conservation equation per se.It's assumed that the incompressible condition is valid including free ice surface points.Second equation stems from the stress-free condition and can be written in terms of stress deviator components.The pressure exclusion leads to the equation at ice surface where Like for the mass conservation equation the mechanical equilibrium equation can be extended to an ice surface.But opposite to the continuity equation the mechanical equilibrium equation is not valid apart from the stress-free condition Eq. ( 6) at ice surface points.So, the extension of the mechanical equilibrium equation to ice surface ∂x . (7) The ISMIP-HOM results have revealed that described approach provides the numerical stability in the finite-difference model at the sparse grids in the experiments A and E (Pattyn et al., 2008).Furthermore, low spread of results for the full Stokes models proves that the approach provides second order of the boundary condition numerical approximation.
The surface temperature was applied as the ice surface boundary condition for the heat transfer Eq. ( 5).This temperature depends on elevation (Makarevich, 1969) as where T s0 (t) was taken equal to −2 • C in the steady state experiments and time dependent annual air temperature in the other experiments.The value −2 • C and the elevation 4160 m above see level in Eq. ( 8) were chosen based on the ice temperature measurements at the Gregoriev Ice Cap summit (Arkhipov, 2004).

Basal boundary conditions
Borehole temperature measurements at the Gregoriev Ice Cap summit in 2004 have shown that the ice temperature particularly doesn't depend on vertical coordinate at the depth larger than 10 m and negative (about −3 • C at the summit).It proves that, firstly, the ice is frozen to the bed (v b =0) and, secondly, the geothermal heat flux is close to zero.increase.But no-slip basal boundary condition in couple with the incompressible equation in the downstream hampers the upstream increase.Vise versa, upstream forces the ice flow in the downstream but the ice flow velocity is relatively small due to no-slip basal boundary condition and is equal to zero at the lowest ice front point.Furthermore, the contrariety leads to instability of the numerical solution of the diagnostic equations Eq. ( 2).It should be underlined that the problem is not purely in the numerical solution but rather in the mathematical statement of the problem.So, it's supposed that the basal drag in the downstream can arise at some environmental changes.The slip boundary condition was applied only at the glacier tongue which was considered as the flowline domain about 600 m upstream from the lowest front point.Basal sliding was introduced through the linear friction law (MacAyeal, 1989(MacAyeal, , 1992)): The basal boundary conditions Eq. ( 9) are the small basal slopes approximation.The common linear friction law can be written as The first equation from Eq. ( 10) follows from the equation σ i k n k =−K f r v i , that is valid only for tangential to the bedrock components (n is the unit normal vector pointing into the bedrock).So, it's incorrect apart from the second equation, which means that the normal to the bedrock basal ice flow velocity is equal to zero.The ISMIP-HOM results have shown significant discrepancy between the different models including full Stokes models, when the friction coefficient suddenly jumps from zero to infinity and vice versa (Pattyn et al., 2008).At this paper is assumed that Introduction

Conclusions References
Tables Figures

Back Close
Full the no-slip/slip transition occurs gradually.From the mathematical point it means the continuity variation from infinity to a constant value of the friction coefficient along the glacier tongue flowline.So, the friction coefficient can be expressed in the form of the Laurent series where X f r is the boundary of no-slip/slip transition.In practice, few main terms of the series were considered.
The basal boundary condition for the heat transfer equation is assigned by geothermal heat flux and basal friction (Pattyn, 2000(Pattyn, , 2003)): 3.3 Boundary conditions at the summit and at the front point The boundary conditions at the summit (x=0) are u=0, σ xz =0, ∂h s ∂x =0.The boundary conditions at the front point (x=L) are v =0 in the case of no-slip boundary condition over the whole flowline domain and the appropriate mass balance equation for the prognostic Eq. (3).In the case of the slip boundary condition we considered the linear interpolation of the depth averaged horizontal velocity from N−2 and N−1 nodes to the last Nth node and the w=0 approximation.

Numerical solution
The problem was solved by the finite-difference method.The solution is based on the two-step algorithm (MacAyeal et al., 1996).First step involves solution of the nonlinear diagnostic equations by an iterative procedure (Hindmarsh and Payne, 1996) for given flowline domain.The second step involves solutions of the heat transfer equation Introduction

Conclusions References
Tables Figures

Back Close
Full and the prognostic equation to obtain the ice temperature and ice thickness distribution, respectively, at the next time step.To avoid the mass conservation problems the staggered-grid convention was used (MacAyeal et al., 1996).

Coordinate transformation
To implement the finite differencing an arbitrary flowline domain is transformed into rectangle by the coordinate transformation.The transformation is performed by replacing the vertical coordinate z to the dimensionless one: ξ=(h s −z)/H.The forms of the equations are modified after the coordinate transformation.In particular the second equation from the system Eq.( 2) in variables x, ξ can be rewritten as where The above mentioned alteration of Eq. ( 13) to ice surface points requires to select the shear stress component σ xz and to replace it with α(2σ xx +σ yy ) in accordance with stress-free surface condition Eq. ( 6).
The terms of Eq. ( 13) are grouped within certain order.The terms of the fist and second lines in Eq. ( 13) represent the longitudinal stress gradients (2 ∂x and ∂σ yy ∂x ) of the mechanical equilibrium equation from the system Eq.(1).They are unchangeable at ice surface points.
The third, fourth and fifth lines of Eq. ( 13) are 2ξ , where ∆ξ is the vertical grid size and upper index "1" means an ice surface grid point.
The integral terms of Eq. ( 13) are equal to zero at ice surface points.The numbers of grid nodes were taken equal to 71 and 31 in the horizontal and vertical directions, respectively.The time step was taken equal to 1 year.

The steady state experiments
The surface mass balance measurements were carried out at the Gregoriev Ice Cap flowline in 1987 and 1988 (Mikhalenko, 1989, Fig. 3).The linear approximation of the averaged values for two years is employed at this paper as the referenced surface mass balance (M s ).Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
The previous investigation of the Gregoriev Ice Cap (Nagornov et al., 2006) was focused on the subsurface stress variation due to seasonal climate changes, but a longterm prognostic experiments to evaluate the enhancement factor in flow law (m, Pattyn, 2000Pattyn, , 2003) ) by comparison the modeled and observed surface profiles hadn't been carried out.Authors took in to account that the stress deviator components are independent from the enhancement factor because if ice flow velocity increases p times then the viscosity (η) decreases p times.The choice of the enhancement factor value was based on the ice surface flow velocity measurements (Vinogradov, 1962).But, the experimental data are very poor and probably high velocity variability close to ice front hadn't been taken into account.The model ice flow velocity is in concordance with the experimental data for relatively small enhancement factor which is evaluated in the range 0.15-0.3.
The results of steady state experiments, when the model runs forward in time until a steady state for the reference surface mass balance is achieved, are represented at Fig. 4.
In the case of small enhancement factor (m=0.15) the glacier front retreats during about 220 years, but ice thickness grows upstream from x≈2.3 km.And later the increased ice mass in the upper part of the glacier return the retreat to advance.Finally, the steady state glacier length is larger than the initial one (3.7 km) and is about 4.14 km.The ice thicknesses along the flowline are 1.5-2 times higher than observed ones at x<2 km (Fig. 4).
On the other hand, the distinction between the modeled steady state and observed surface profiles is noticeably less for m=1, while ice thicknesses at x>2.5 km smaller than the measured values (Fig. 4).Some disagreements in the profiles at x<2 km can be explained by insufficient account of the lateral ice flow in the model.negative at elevations from about 4600 m a.s.l.(Fig. 3).The full deglaciation time for the displacement ∆M s =−0.6 m/yr is about 1200 years (Fig. 6).

Non-steady experiments
In the non-steady experiments the periodical climate variability was considered, when the annual air temperature varies as where t p is the climate period, ∆T a is the air temperature amplitude which was taken equal to 0.3 It was supposed that the annual air temperature and the ice surface mass balance vary proportionally and the dependence can be empirically expressed as where M s0 is the reference mass balance.
The mass balance and air temperature amplitude relation is suggested unknown, but the amplitudes have the opposite signs.It means that warming/cooling decreases/increases ice surface mass balance.
The experiments were carried out for different mass balance amplitudes (∆M s ) in the range −0.6−0.2 m/yr at the climate period t p =500 yr (Fig. 7) and for different climate periodicity at the fixed amplitude ∆M s =−0.4 m/yr (Fig. 8).Introduction

Conclusions References
Tables Figures

Back Close
Full Screen / Esc Printer-friendly Version

Interactive Discussion
The basal sliding arises at certain time intervals when the ice mass grows.We considered the simple criteria of the no-slip/slip transition in the glacier tongue base.It was formulated for the basal shear stress as max X tr <x<L (σ xz ) b >σ cr , where σ cr was considered in the range 1-1.2×10 5 Pa.The points of the transitions at the curves look like the derivative break points (Figs. 7 and 12).Really, the transitions are more complex and probably occur smoothly.The friction coefficient values were considered timeindependent.More intensive ice mass growing requires more rapid glacier advance and, hence, smaller friction coefficient values.The values were chosen accordingly Eq. ( 11) and closely to a maximum as it possible to avoid the above mentioned instability in the diagnostic system solution.So, in the experiments the friction coefficient depends on the mass balance amplitude ∆M s (Fig. 9).The glacier length oscillations periods are in concordance with climate periods.Despite the diagnostic problem nonlinearity the glacier length amplitude versus the mass balance amplitude is close to the linear relationship (Fig. 10).While the deviation for ∆M s =−0.6 can be explained by the bedrock topography inhibition (Fig. 11).
The tree-rings data collected at the Tien-Shan (ISTC, Project No. 2947) allowed to define the short-term climate periodicity in the region.The 12, 60, 110 yr periodicity had been defined and with assumption of existence of long-term periodicity, which was taken equal to 800 yr, the superposition of the equal amplitudes (∆T a =0.2 • C and ∆M s =−0.2 m/yr) leads to the glacier length history as shown at Fig. 12.According to the glacier length amplitude versus mass balance periodicity (Fig. 8) short-term climate periodicity insignificantly affects the glacier length changes which are governed by the long-term component.
It can be suggested, that the glacier can advance by the snow accumulation at the surrounding ground when the surface mass balance increases.In the considered experiments the mass balance at the front point (x=L) is negative with the exception of some time intervals in the cases of |∆M s |>0.4 m/yr (Fig. 13).From the point, when the mass balance at x=L becomes positive, to the point, when the intensive glacier advance begins due to sliding, the firn layer with about 5 m thickness can accumulate

Conclusions References
Tables Figures

Back Close
Full from x=L to the distance about 700 m (∆M s =−0.6 m/yr).The firn thickness was estimated in the assumption that the mass balance is the same at ice surface and ground.But, the albedo is different and, it seems, firn thickness should be smaller.So, we supposed that the firn layer can't sufficiently affect the basal shear stress at the glacier tongue at the point of the basal sliding beginning.
The melting rate at the base can be estimated as 5 Pa and u b ≈10 m/yr which were observed in the experiments.So the ice surface mass balance dominates in the prognostic equation.
As the sliding occurs only at the glacier tongue base, the glacier length deviations are insignificant (<1%) for considered friction law approximations by Eq. ( 9) and by Eq. ( 10), respectively (Fig. 13).For comparison, more significant distinction in the glacier length changes is observed for the different vertical grid sizes (N ξ =11 and N ξ =31), when the relative deviation achieves 2.5%.The vertical resistive stress effect accumulates in time, but also leads to insignificant length deviation (<2%) in the case of periodically climate changes.

The relationship between glacier length and annual air temperature
We investigated the correlations between glacier length and annual air temperature changes within the dependence was suggested by Oerlemans: Basing on the non-steady experiments it's clear that there are no unique parameters α, β and γ in Eq. ( 16), because the glacier length amplitude significantly depends on the climate periodicity at fixed air temperature amplitude and mass balance amplitude (Fig. 8).But it seems the main (long-term) climate period can be successfully derived from glacier retreat/advance history based on Eq. ( 16) and taking into account the dependences of α, β, γ versus the periodicity.

TCD Introduction Conclusions References
Tables Figures

Back Close
Full For given temperature history and glacier length changes obtained by ice flow modeling the parameters α, β, γ satisfy the minimum of discrepancy between left and right members of Eq. ( 16).The discrepancy minimum conditions lead to the linear algebraic equations about α, β, γ: The major drawback of the air temperature reconstruction from glacier length history according to Eq. ( 17) is the time derivation ( d L d t ) sensitivity with respect to (i) bedrock undulation, (ii) features of glacier advance (with basal sliding or without it), (iii) front shape reorganization prior to retreat/advance beginning.We tried to overcome the problem by smoothing procedure which was applied to the glacier length history (Fig. 15).At the first glance, the above mentioned features of the glacier flow are invisible in the glacier length history, but the success of the temperature reconstruction significantly depends on the "degree" of the glacier length history smoothing (Fig. 16).
On the other hand, the smoothing can deprive the possibility of the reconstruction of short-term temperature changes which weakly appear in the glacier history (Figs. 12     and 19).So, it's difficult to point optimal "degree" of glacier length history smoothing and the temperature paleoreconstruction requires coupling of independent methods such as moraine sediments analysis and, for example, the δ 18 O data analysis or glacier surface temperature paleoreconstruction based on the ice borehole temperature measurements.
It seems, another way to overcome the problem is to introduce the initial phase ϕ 0 into temperature-time relationship, i.e. to replace Eq. ( 16) with a new one in the form T a (ωt + ϕ 0 ) = α + βL(t).

Conclusions References
Tables Figures

Back Close
Full The γ versus periodicity (Fig. 18) shows that the initial phase ϕ 0 sufficiently depends on the temperature harmonic oscillations frequency.The initial phase dispersion supposes that the reconstructed temperature can be a priori represented in the form of a part of harmonic series.So, such approach application seems more poor beside the one based on the Eq. ( 16).
The dependences of α, β, γ versus mass balance amplitude (Fig. 17) can be derived by ice flow modeling and the data can be collected for the future climate paleoreconstructions based on moraine sediments analysis in the different regions of the Earth.

Conclusions
2-D flowline ice cap modeling results show that deglaciation of the plane tops of Terskey Ala-Tau is observed for the significantly smaller values of the ice surface mass balance than the reference mass balance based on the experimental data which were obtained in 1987 and 1988.In particular, mass balance at the Gregoriev Ice Cap surface should be smaller by about 0.6 m/yr than the reference one, and for the reference mass balance the glacier will be almost unchanged.
Hence, the correct forecast of future glaciation/deglaciation of the plane tops requires the confirmation of the reference mass balance and resumption of monitoring of the ice surface mass balance changes.
2-D flowline ice cap modeling results show that, in general, annual air temperature and glacier length history correlations for the harmonic climate changes are well described within linear dependence of annual air temperature versus glacier length and time derivation of the length.But, in particular, the dependence is sensitive to nature of the ice flow as well as to bedrock undulations through the time derivative of glacier length.So, the air temperature history reconstruction requires the smoothing of the glacier length history.On the other hand, the smoothing can lead to losses of data of short-term climate changes.The problems indicate that correct climate reconstruction requires integration of independent methods such as moraine sediments analysis, specific latent heat of fusion 3.35×10 after the pressure exclusion.The manipulations lead to two equations in the 3-D model or one equation in the 2-D flowline model, respectively, in terms of stress deviator components.
It's suggested, that the surface mass balance can increase significantly and steep in time as a result of the climatic variability.Such climatic scenario leads to the ice stream Screen / EscPrinter-friendly Version Interactive Discussion Screen / EscPrinter-friendly Version Interactive Discussion So, basing on the steady state experiments, there are no reasons to correct the flow-law rate factor values(Paterson, 1994)  for the plane tops glaciers.Visible glacier retreat/advance didn't occur for the referenced surface mass balance (m=1).Hence, deglaciation of the plane tops should arise at the smaller values of the decreased the referenced M s by the step −0.1 m/yr until the full glacier disappearance happened.The steady state profiles have been achieved for the parallel displacement (∆M s ) of the referenced mass balance by the values −0.1-−0.5 m/yr, and are represented at the Fig.5.The full deglaciation occurs if the mass balance is less than the referenced one by −0.5 m/yr, then the mass balance becomes

Fig. 5 .Figure 6 .
Fig. 5.The ice cap steady state surfaces for different ice surface mass balances.

Fig. 6 .
Fig. 6.The ice cap length changes in the steady state experiments.

Fig. 13 .Figure 16 .
Fig. 13.Mass balance history at the front point in the cases of different mass balance amplitudes.

Fig. 16 .
Fig. 16.The annual air temperatures derived from the ice cap length history.

Table 1 .
analysis and glacier surface temperature paleoreconstruction based on the ice borehole temperature measurements.The ice flow modeling allow to obtain the parameters of the temperature-length dependence for different climate histories and the data can be collected for the future forecasts or climate reconstructions.List of symbols.