Influence of anisotropy on velocity and age distribution at Scharffenbergbotnen blue ice area

We use a full-Stokes thermo-mechanically coupled ice-flow model to study the dynamics of the glacier inside Scharffenbergbotnen valley, Dronning Maud Land, Antarctica. The domain encompasses a high accumulation rate region and, downstream, a sublimation-dominated bare ice ablation area. The ablation ice area is notable for having old ice at its surface since the vertical velocity is upwards, and horizontal velocities are almost stagnant there. We compare the model simulation with field observations of velocities and the age distribution of the surface ice. No satisfactory match using an isotropic flow law could be found because of too high vertical velocities and much too high horizontal ones in simulations despite varying enhancement factor, geothermal heat flux and surface temperatures over large ranges. However, the existence of a pronounced ice fabric may explain the observed present-day surface velocity and mass balance distribution in the inner Scharffenbergbotnen blue ice area. Near absence of data on the temporal evolution of Scharffenbergbotnen since the Late Glacial Maximum necessitates exploration of the impact of anisotropy using prescribed ice fabrics: isotropic, single maximum, and linear variation with depth, in both two-dimensional and threedimensional flow models. The realistic velocity field simulated with a noncollinear orthotropic flow law, however, produced surface ages in significant disagreement with the few reliable age measurements and suggests that the age field is not in a steady state and that the present distribution is a result of a flow reorganization at about 15 000 yr BP. In order to fully understand the surface age distribution, a transient simulation starting from the Late Glacial Maximum including the correct initial conditions for geometry, age, fabric and temperature distribution would be needed. This is the first time that the importance of anisotropy has been demonstrated in the ice dynamics of a blue ice area and demonstrates the need to understand ice flow in order to better interpret archives of ancient ice for paleoclimate research.


Introduction
Blue ice areas (BIAs) make up about 1 % of the surface area of Antarctica (Bintanja, 1999).They are usually formed where removal of snow fall by wind exceeds precipitation (Crary and Wilson, 1961), though low-lying regions allow surface melt to be the dominant process.The Scharffenbergbotnen BIA (see Fig. 1), as with most BIAs, is located in a region where wind speeds are higher than normal for Antarctica.High winds are generally because of steep gradients promoting katabatic flows, or in mountain areas where local wind fields determined by the geometry of mountains and ice surface lead to accelerated flow (Takahashi et al., 1992).Ice flow then brings deeper ice to the surface, exposing aged solid blue ice (Naruse and Hashimoto, 1982;Whillans and Cassidy, 1983;Azuma et al., 1985;Bintanja, 1999).The lower albedo of bare ice surfaces relative to snow-covered areas enhances ablation in BIAs (Bintanja and Reijmer, 2001).
Published by Copernicus Publications on behalf of the European Geosciences Union.
Consequently, isochrones are inclined relative to the surface -from sub-horizontal near the equilibrium line to near vertical at the stagnant ends of the valley (Naruse and Hashimoto, 1982;Whillans and Cassidy, 1983).This makes BIAs an ideal location to obtain "horizontal ice cores" -that is, collecting samples for analysis along a shallow surface trench rather than traditional deep drilling -which may contain a record of climate record spanning the Holocene and perhaps longer (Moore et al., 2006;Sinisalo and Moore, 2010).Presently, however, the problems in dating and interpreting climate records from horizontal ice cores have largely precluded their usefulness as climate archives (Sinisalo and Moore, 2010).Chief among these difficulties is an understanding of the long-term stability of the blue ice as the surrounding ice sheet changes through a glacial cycle.
Ice sheet elevation changes at the glacial termination are likely to have been most pronounced in the nunatak areas, such as Scharffenbergbotnen, a few hundred kilometers from the coast (Pattyn and Decleir, 1998), since this is the transitional area between cold-based and warm-based ice flow.The surface of the major part of the East Antarctic Plateau ice sheet may have been about 100 m lower in the last glacial than at present (Jouzel et al., 1989;Pattyn, 1999;Ritz et al., 2001).In contrast, the surface elevation at the margins of the plateau may have been hundreds of meters higher (Näslund et al., 2000;Hättestrand and Johansen, 2005).
Understanding the present-day ice dynamics of BIAs is challenging.Several attempts have been made using simplified flow models constrained by limited data (Naruse and Hashimoto, 1982;Azuma et al., 1985;Van Roijen, 1996;Grinsted et al., 2003;Sinisalo et al., 2004Sinisalo et al., , 2007;;Moore et al., 2006).However the complex geometry, usually including large ice thickness variations along the blue ice field, suggests that a three-dimensional, higher order or full-Stokes model is needed to produce realistic ice flow simulations.Scharffenbergbotnen is the best-studied BIA with a long history of observations including ice depth, surface velocities, mass balance, ice dating, ice temperature measurements, and both vertical and horizontal ice core archives of ice chemistry and water isotopes (the data sets are summarized in Sinisalo and Moore, 2010).Hence this BIA is a suitable test bed for numerical simulation using advanced flow modeling, which could provide insights on the dynamical evolution of the region since the last glacial.
After presenting the input data in Sect.2, we introduce our model and the chosen initial and boundary conditions (Sect.3).The simulation setup as well as the basic concepts of the applied anisotropic rheology follow in Sect. 4. The different approaches in terms of simulations are presented in Sect. 5 and then discussed in detail in Sect.6, from which the conclusions (Sect.7) with respect to the dynamics and age distribution at Scharffenbergbotnen are deduced.
2 Study area: Heimefrontfjella, Scharffenbergbotnen Scharffenbergbotnen (SBB), northwest Sivorgfjella (74 • S, 11 • W; 1200 m a.s.l.), is a closed, glaciated valley with inflow from the surrounding ice sheet.The head of the valley acts as a dam on the main flow.There are two separate BIAs (light-gray in Fig. 1) in the valley, however this work focuses on the innermost one.Today, ice flows into the valley mainly from the wide north-northwestern entrance, generally with an eastward-directed flow (see stake 6 in Fig. 1) and, to a very limited extent, over the mountains in an ice fall at the eastern end of the valley (close to stake 19 in Fig. 1).Geomorphological evidence (Hättestrand and Johansen, 2005) suggests that during the Late Glacial Maximum (LGM) the ice surface in the valley was 200-250 m higher than today, while the ice sheet outside the valley (i.e., outside the area confined by the nunataks) was only 50-150 m higher.The LGM ice sheet was able to overflow some of the passes along the valley walls and ice entered the valley from several more locations than today.As the elevation of the surrounding ice sheet decreased after the LGM (e.g., Hättestrand and Johansen, 2005), the ice overflow became insignificant and the ice flow inside the valley decoupled from the main ice sheet flow.
A geological map of the Scharffenbergbotnen area is available from aerial photogrammetric surveys (1993, scale 1 : 25 000) as is a digital elevation model of the surface (100 m resolution) based on airborne surveys in 1985/86 over the whole Heimefrontfjella area (both copyright of Bundesamt für Kartographie und Geodäsie, Frankfurt am Main, http://www.bkg.bund.de).Ice thickness data (100 m resolution) comes from a radio-echo sounding campaign in 1987/88 and generally increases from the flanks of the nunataks surrounding the valley towards the center of the valley but with a maximum thickness of 1050 m at the northwestern entrance (Herzfeld and Holmlund, 1990).
Several ice cores have been drilled in the area that provide data on the age of the near-surface ice.A 52 m long vertical core was drilled in summer 1997/98 in the inner part of the valley (V in Fig. 1).The uppermost 45 m section has been dated to 10 500 (+700/−300) calendar years BP (van der Kemp et al., 2002), i.e., age calibrated for atmospheric 14 C variations.A 100 m horizontal ice core (the surface of the ice was sampled along a trench) was collected further upstream in 2003/04 (at H in Fig. 2).The middle part of this ice core was dated to 4426 ± 215 calendar years BP (Sinisalo and Moore, 2010).A longer horizontal core (2.7 km) was collected in 2006/07 (the complete blue line in Fig. 2).
The age of the ice has been estimated at several locations (Table 1).It should be pointed out that there are some uncertainties in age observations, especially of ages in the valley.Some of the older measurements were made with a technique based on the 14 C age of gas trapped in air, and this is known to be subject to error caused by in situ production of 14 C by cosmic rays that penetrate surface ice layers (van der Kemp et al., 2002).Newer data come from measurements on  .The inset on the upper left side is a close-up of the area around the horizontal ice core.Age measurements are shown as pentagons with the coloring according to the age legend with a distinction between reliable data (filled pentagons with black outlines) and data points with larger errors (empty pentagons with outlines in the color corresponding to their value) -see Table 1.   1, Fig. 2).Other locations were sampled, including the relatively darker band (orange feature in Fig. 2) in the ice at position b15, but no other samples contained enough carbon to produce a reliable date.Van Roijen (1996) developed a method for dating blue ice by measuring 14 C concentration in air trapped in the ice.In Van Roijen's method, 14 C depth profiles in blue ice are translated into carbon ages with a correction made for 14 C produced in situ, requiring an ice core to be analyzed to about 50 m depth, and that results in the age at point V, Table 1.Additional measurements are available from 14 C analysis in 10 other locations; see Fig. 2 (numbered empty pentagons).These 14 C ages had large uncertainties of up to several thousands of years (Van Roijen, 1996) since they were collected from shallow ice cores without correction of in situ 14 C production (van der Kemp et al., 2002).The results of the 14 C measurements from the shallow ice cores (Van Roijen, 1996) have been corrected for atmospheric 14 C variations and calibrated to calendar ages using OxCal, an online calibration tool (http://c14.arch.ox.ac.uk).
The errors given in Table 1 correspond to the standard deviations given by the calibration.Since these measurements are much less reliable than the two measurements obtained at V and H, they are only used for approximate age estimation.In Sect.5.3 we introduce a flow-line model to reduce the full three-dimensional solution complexity and to study principal influences from anisotropic rheology.The flow line was designed to coincide with the longer horizontal ice core (see Fig. 2).Inspection of aerial photos and maps, together with in situ observations of the blue ice led to identification of a possible flow line on the valley.The main features used to identify this flow line were the visible bands on the ice surface (see e.g., Sinisalo and Moore (2010) -and several are also visible on the geological orthophoto map).We selected the route so that it passed through the apex of these curving structures as we supposed the ice velocity was fastest along those.The flow line terminated in the inner valley at a broad (3 m wide) band of dark-colored ice that was about 1 m above the surrounding ice, and which we interpret as the confluence of the ice flow from the main entrance to the valley and the much smaller inflow from the ice fall.This convergent flow maintains a higher ice elevation on the band despite being subject to enhanced ablation due to its low albedo.The flow line defined the sampling location of the horizontal ice core.Note that this flow line is different from our previously suggested flow line (Grinsted et al., 2003), when we believed that flow originated mainly outside the valley for the inner blue ice region.Now we think that the main source of ice for the northwest part of the inner BIA is the snow accumulation region between the two BIAs -east of stake 4 north of stake 22 (see Fig. 1) -which has an especially high accumulation rate region near the northeast valley wall.
A geothermal heat flux of q geo = 60 mW m −2 was estimated at the position of SBB based on Shapiro and Ritzwoller (2004) -extracted from the SeaRISE project (http:// websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica).Present-day mean annual temperatures on the blue ice area are estimated as −20 • C based on observations from IMAU AWS 7, which operated from December 1997to January 2003(C. H. Reijmer, personal communication, 2010).Since the LGM surface temperatures have varied by about −8 to +3 K (in the Holocene climate optimum) relative to those at present (Sato and Greve, 2012).Surface elevation changes of 250 m would modify these by only a degree or so, and this range of variability is consistent with stable water isotope variability on the BIA (Sinisalo et al., 2007).2001/02, 2002/03, 2003/04 and 2006/07.Not all stakes were resurveyed, and some could not be relocated at all.All stakes were surveyed by differential GPS using a base station tied to a reference marker located on bedrock near the Svea station (Fig. 1).Baselines of less than 20 km ensured centimeterrelative precision.The largest sources of error were variable stake lean or measurement during windy conditions; however this error was generally much less than the motion of the stakes between surveys over different seasons.Stakes were lost on the high accumulation snow hill between the two BIAs (Fig. 1).For each stake the velocity calculated over the longest available period has been used in our study after checking that the measurements over shorter periods are consistent with the long-term trend.Some results published in earlier articles have been improved upon here.Horizontal velocities are very low; the ice is hardly moving in this area compared with ice flow in the surrounding ice sheet.Over large areas of the BIA surface velocities are below 1 m yr −1 .The stake velocities are shown in Fig. 1 and Table 2.In addition we used some earlier measurements (Ber. Polarforsch. 86, 1991) surveyed using theodolites (see Fig. 1) within one single season that are less reliable than GPS surveys over longer periods.

Model
By the nature of its high viscosity, the flow of ice is described by the Stokes equation, which neglects inertia forces and sets the specific driving force, namely density times gravity, ρg, in balance with the divergence of the Cauchy stress, divσ , divσ + ρg = 0. (1) The Cauchy stress tensor, σ = τ −pI , can be divided into the deviatoric stress, τ , and the isotropic pressure, p = −trσ /3.The standard way of treating ice in glaciers and ice sheets is to assume incompressibility, which is justified if the compression of firn has a negligible influence on the global dynamics of the system.Although firn compaction in the high-accumulation area between the two BIAs could have an influence on the local dynamics, in the case of blue ice, where by definition we have no firn layer to account for, this certainly holds and the conservation of mass reduces to divu = trε = 0. (2) The strain rate tensor ε can be deduced from the velocity vector, u = (u, v, w) T , by If treated as an isotropic material, ice rheology is given by a Norton-Hoff power law, with the power law exponent n = 3 in Eq. ( 5) known as Glen's law in glaciology, which collinearily links the deviatoric stress τ with the strain rate tensor ˙ : where the effective viscosity η is defined as In Eq. ( 5), ε2 e = tr(ε 2 )/2 is the square of the second invariant of the strain-rate tensor.Under applied shear, the ice grains rearrange their c axes from a random anisotropic distribution into a more energetically favorable type of fabric (Duval, 2000) and hence the rheology will deviate from isotropic behavior.To a limited extent, the enhancement factor E can be used to mimic anisotropy effects (e.g., Ma et al., 2010).A = A(T ) is a rheological parameter that depends on the ice temperature T relative to the temperature melting point via an Arrhenius law (Paterson, 1994).
The temperature, T , is obtained from the general balance equation of internal energy and reads where κ = κ(T ) and c v = c v (T ) are the heat conductivity and specific heat of ice, respectively.The last term in the heat transfer equation represents the amount of energy produced by the viscous deformation, which because of the low shear rates in our application has no significant contribution and hence is neglected.
The complete set of equations is solved using the finiteelement method (FEM) using the open source software package Elmer/Ice (Gagliardini et al., 2013; http://elmerice.elmerfem.org),which has already been applied to glacier simulations of different kinds (Zwinger et al., 2007;Zwinger and Moore, 2009;Zhao et al., 2014), and we use similar numerical settings as discussed in Martín et. al. (2014).
The age of the ice is governed by the advection equation Because of its purely advective nature, we chose to solve Eq. ( 7) using a semi-Lagrangian scheme integrated over a period of up to 15 000 yr using the steady-state velocity field obtained with a thermo-mechanically coupled solution.Exactly the same semi-Lagrangian method was first applied and explained by Martín et. al. (2014).In Sect. 5 we chose to present results taken after 15 000 yr, as this corresponds to the time span since the middle of the last deglaciation in Antarctica (EPICA community members, 2004).

Anisotropic rheology
Until recently, anisotropy in Elmer/Ice was restricted to the linear general orthotropic flow law  (GOLF; Gillet-Chaulet et al., 2006) or a colinear CAFFE model based on an anisotropic flow enhancement factor (Seddik et al., 2008(Seddik et al., , 2011)).In more recent work (Martín et al., 2009;Ma et al., 2010), the orthotropic law has been extended to a nonlinear form by adding an invariant into the anisotropic linear law.As we will show later, one of our main findings was the necessity to introduce anisotropic flow behavior to understand the ice dynamics of the Scharffenbergbotnen valley and match it qualitatively with the combined data of observed velocities, surface mass balance and surface age distributions.We did this by applying the model introduced and discussed in detail by Martín et. al. (2014).Then the noncovariant formulation of anisotropic flow Eq. ( 4) becomes a general tensor relation, which in index notation can be written as The 36 independent components of the viscosity tensor η are functions of the polycrystalline fabric.Mathematically, this dependency is expressed in terms of the second-and fourthorder orientation tensors a (2) and a (4) , respectively, defined as with c being the c-axis unit vector and <> the volume average.In order to express the fourth-order tensor in terms of the second-order tensor (Advani and Tucker, 1990), a closure relation is provided (Chung and Kwon, 2002;Gillet-Chaulet et al., 2006).

Boundary and initial conditions
The governing equations have to be accompanied by boundary conditions for the field variables.All our simulations are performed on a fixed geometry, reducing the otherwise transient kinematic free-surface condition to The Cryosphere, 8, 607-621, 2014 www.the-cryosphere.net/8/607/2014/based on steady-state assumptions.In this way the surface net accumulation, a s , becomes part of the solution of the system rather than an additional boundary condition to it.By neglecting atmospheric pressure gradients and shear forces exerted by the atmosphere, the dynamic boundary condition at the free surface reduces to a vanishing Cauchy stress vector, where n is the outward pointing surface normal, The remaining variables to be set at the free surface are the present annual mean temperature -taken to be T z s = −20 • C -as well as the age of the downwards advected ice Because of the special hyperbolic nature of Eq. ( 7), this condition is set only if u • n < 0. In other words, only if the velocity is pointing into the glacier do we set Eq. ( 14).
Even with the relatively large heat flux of q geo = 60 mW m −2 we used, only a small area of the bedrock, beneath the deepest ice at the very outer end of the valley (where it has a minor impact on the inner BIA dynamics), reaches pressure melting point.This allows for setting a fixed boundary condition, u = 0, for the velocity at the bedrock.This, in consequence, prohibits a steady-state solution of Eq. ( 7), as it would lead to infinite age, → ∞, at the bedrock.As an initial condition for the age distribution, (x, t = 0), we assumed plug flow and applied the analytic solution of Lliboutry (1979) using a constant artificial melt-rate of 1.0 × 10 −6 m yr −1 to allow for steady-state conditions to be applied.

Simulation setup
The velocity field resulting from the Stokes equation is an instantaneous response to a given viscosity distribution and not directly influenced by temporal changes of the geometry.On the other hand, temperature and fabric evolution as well as the age-depth distribution depend on the flow history of the glacier and would demand a transient simulation starting from the LGM and extending to the present day.In particular, by our simulations it turned out that the age-depth distribution -even for fixed geometry -is not in a steady state and by the purely convective nature of Eq. ( 7) demands a transient approach even starting from before LGM to get the correct initial age-depth distribution.Such a simulation demands -the LGM geometry (initial condition for z s ); -the LGM age and fabric distribution, LGM and a (2) ij,LGM (or spin-up simulation starting one glacier cycle earlier); -the accumulation-ablation pattern during the whole simulation (a s (x, t) in Eq. 10); -the surface temperature history from the LGM to present T s (x s , t).
As only limited geomorphological information on LGM geometry is available (Näslund et al., 2000;Hättestrand and Johansen, 2005), at least the first two of these items are largely undetermined.Additionally the computational effort required to run a full-Stokes model for such a long period, with relatively unconstrained geometrical evolution, renders a prognostic simulation from LGM a highly unrealistic approach.Consequently, we focused our efforts on studying the effects of anisotropy on ice dynamics within the fixed present-day geometry assuming a steady thermal state and prescribed fabrics.
For all simulations we assume steady-state ice-flow conditions and calculate the evolution of the age distribution from a typical ice-sheet configuration to a blue-ice area.To that end, we calculate the steady-state temperature and flow distribution under the present conditions.Using the velocity solution from the steady-state run, the age Eq. ( 7) is integrated from an initial vertically layered plug-flow distribution using a transient run, which -motivated by the values of the most reliable age measurements at the surface -starts at 15 kyr BP.
Based on the data sets presented in Sect.2, a twodimensional, unstructured footprint mesh, covering the glaciated area of SBB in the horizontal plane was created using the open-source meshing software Gmsh (Geuzaine and Remacle, 2009).Horizontal resolution is 30 m in our main region of interest (close to the ice cores), 50 m in the rest of the valley and up to 100 m closer to the boundaries.We extruded this mesh in the vertical direction in 13 levels.The resulting mesh has about 61 000 nodes that are connected in 54 000 linear elements.
Elmer/Ice uses the finite-element method to discretize the governing equations.The resulting system matrices for the solution of the temperature field as well as the Stokes equation (velocities and pressure) are pre-conditioned with an incomplete LU decomposition and thereafter solved using the generalized conjugate residual (GCR) method.Because of the nonlinearities introduced by the material parameters (heat conductivity and capacity as well as viscosity) each solver for itself has to iteratively solve a linearized system, obtained by a fixed-point iteration for resolving these nonlinearities.Due to their mutual coupling we also have to iteratively solve for temperature and the Stokes equation to obtain a converged solution for the steady state.
Three different setups were tested to study the principal influence of anisotropic rheology on the flow and the agedepth distribution: (i) Isotropic flow using Glen's flow law.
(ii) A depth-dependent, piecewise linear anisotropic fabric distribution starting from isotropic and developing with the local ice thickness h(x, y) and the distance from the free surface, d.
(iii) An over-the-depth fixed single maximum anisotropic fabric: Note that Eq. ( 15) is defined for three dimensions, and also applied for the flow-line model.The isotropic and the single maximum case define the extreme configurations of the fabric that may occur and comparison with measured values will show how pronounced the orientation of the c axis is to be expected.For all 3-D runs we further compute the corresponding surface accumulation distribution, a s , by inserting the diagnostically computed velocity field (u s , v s , w s ) T as well as the given surface geometry, z s into the right-hand side of Eq. ( 11).Results of this comparison are listed in Table 3, from which it is clear that the reconstructed computed steady-state surface mass balance obtained with isotropic rheology is off by multiple orders of magnitudes in some places.

Results
This section presents the numerical simulations we conducted.First, we present the isotropic, three-dimensional (3-D) case, which leads to the problems of having too high horizontal velocities in combination with far too old ice at the measured ice core positions (see Sect. 2).We attempted to adjust the ice viscosity to produce a reasonable match to observations by globally tuning the enhancement factor, E, in relation ( 5) as well as changing the internal temperature distribution -and thereby the rate factor, A(T ), in Eq. ( 5).We did this by varying the surface temperature using the constraints given by climate history within realistic values and the geothermal heat flux.However, it was clear that we could not resolve the velocity and age issue.Figure 3 summarizes the characteristics of the flow field computed for different altered parameters in the isotropic flow law (enhancement factor, temperature and mechanical boundary conditions).To that end, we sampled surface velocities of the inner BIA (see Fig. 1), took the average for the absolute values of horizontal and vertical velocity as well as the horizontal components and normed them by the corresponding values of the reference run (isotropic rheology with enhancement factor 1). From this it is clear that the related changes obtained with variations of enhancement factors or boundary conditions in The other column sets represent results from simulations with a single parameter varying from this standard configuration as denoted by labels beneath the columns.BC represents a simulation with no inflow, u = 0, on the otherwise free inflow boundary of the domain.The rightmost set contains the results for the single maximum fabric simulation, and the inset emphasizes differences in horizontal and vertical components between this case and the isotropic case with the reduced enhancement factor (E = 0.1).
an isotropic setup are only able to alter the magnitude of all components simultaneously.Additionally, only a strongly reduced enhancement factor of E = 0.1 would lead to horizontal velocities that are sufficiently small to match the measured order of magnitude.From Fig. 4 it is clear that there has to be a differentiation for the horizontal components of the velocity in order to match the directions.Furthermore, the discrepancy of the computed steady-state surface mass balance with respect to the measured (see Table 3) for the isotropic as well as the depth-dependent case indicate that the vertical surface velocity component also has to adapt in a nonisotropic way.The most plausible remaining explanation for the mismatch of simulated and observed velocities is the existence of a developed ice fabric that alters resistance to flow differently in horizontal and vertical direction.

Simulations using isotropic rheology in 3-D
The first attempt to model the three-dimensional flow at SBB was to deploy the standard isotropic Glen's flow law (Paterson, 1994) 2 for values).Discrepancies in both direction and magnitude (with observed velocities being almost an order of magnitude smaller) can be easily seen.
strongly influence the steady-state velocity field (Zwinger and Moore, 2009), a short prognostic run allowing for the free surface to adjust for these errors was made.Since no major adjustments of the free surface were observed we chose to use the unaltered geometry for the thermo-mechanically coupled steady-state simulations.We use the resulting velocity field and the semi-Lagrangian method to integrate the age advection Eq. ( 7) over a time span of 15 kyr.
The result of this run is presented in Figs. 4 and 5. Comparing this with the measured velocities (Fig. 1 and Table 3) it is immediately clear that -although the regions of ablation (over the BIA) show a relatively good match -simulated values of horizontal surface velocities exceed measured ones by almost one order of magnitude.They also deviate considerably in direction in some places, in particular in the upper part of the inner BIA (stakes S3, S4, 13 and 26 in Fig. 1).Simultaneously, the simulated age distribution does not match the measured dates (see Table 1 and Fig. 2).They neither match the absolute values -they are far too old -nor in terms of the relative distribution; the gradient between the positions of the two most reliable measured values is much too large.

Simulations using prescribed fabrics in 3-D
In order to investigate the possible influence of anisotropic fabric on the flow inside SBB valley, we altered the threedimensional setup of SBB prescribing the anisotropic distributions individually for Eqs. ( 15) and ( 16).The results are depicted in Figs. 6, 7 and 8.
The surface velocities in case (ii) with depth-dependent anisotropic fabric (Fig. 7a) are insignificantly different from the isotropic case (see Fig. 4).The issue of too high velocities (by an order of magnitude) with ill-matching directions remains.Also, the surface age distribution (Fig. 8a) resembles the one obtained with the isotropic flow law (see Fig. 5).
A clear change is observed with case (iii), the single maximum anisotropic fabric, as can be seen in Fig. 3.The velocities (Fig. 7b) are greatly reduced with respect to the isotropic case and now match the directions, with absolute values slightly lower than observed magnitudes.Deeper in the ice, the isochrones become more horizontally aligned  (Fig. 6c), shortening the distance an ice particle has to travel from the accumulation area to the inner valley.For the inner BIA this results in a shifting of the old ice towards the position of the measurement point V (green pentagon in Fig. 8) and hence -at least qualitatively -improving the match with the observed age distribution.The slightly too slow velocities lead to almost stagnant ice down the valley (stake Bi1 in Fig. 1), resulting in too old ice ages there (Fig. 8b).A similar trend can be observed for the computed steady-state surface mass balance (SMB).Whereas the isotropic as well as the depth-dependent fabric lead to negative values of the computed SMB that exceed the measured results by an order of magnitude, the SMB obtained with the single maximum fabric is reduced and now generally matches the observations within their error margins.
The change between the isotropic/depth-dependent anisotropic and the single maximum anisotropic case can be even more clearly seen in the corresponding 3-D stream lines.To that end we ran a stream-line solver (Runge-Kutta) provided within ParaView (Ahrens et al., 2005) backwards from a line at the surface that connects the sample points H and V (see Fig. 1 and Table 1) in order to visualize the origin of ice particles that reach the surface in that region.The results of these post-processing steps are shown in Fig. 9.For the isotropic and depth-dependent anisotropic fabric cases (see Fig. 9a), the origin of most particles is either from the region between the two BIAs or the southern flank on the other side of the surface moraine in the inner part of the valley (depicted in Fig. 1), which in itself is already a very unrealistic result (as stream lines are not expected to cross surface moraines).This picture changes for the single maximum anisotropic fabric distribution (see Fig. 9b), where all stream lines are shifted to originate from areas close to the bedrock from within the outer part of the valley, more resembling the flow line inferred from observations of banding in the valley and the flow line marked in Fig. 1 used for the previous studies (Grinsted et al., 2003).Another clear distinction between the two configurations is the integration time, i.e., the travel time of a particle as a result from the Runge-Kutta method to advect a particle from the surface along the stream line back to that particular position.Naturally, given the faster velocities, these travel times are much faster for the isotropic/depth-dependent anisotropic than the single maximum anisotropic fabric, with the latter quickly exceeding 50 000 yr once reaching the outer part of the valley.

Comparison of 3-D results with 2-D flow-line model
In addition to the complete three-dimensional setup, we also investigated the dynamics using a flow-line model, as there have been earlier attempts to compute the age distribution with such models (Grinsted et al., 2003).Initially chosen to reduce the complexity and the size of the problem, the twodimensional (2-D) flow-line model illustrates the necessity of using a full three-dimensional flow simulation.
Comparing Fig. 10 and Fig. 6a and c shows that the results differ quantitatively for both extreme scenarios, the isotropic and the single vertical maximum fabric.The upwelling region of old ice around 5.5 km along the flow line observed in the 3-D isotropic simulations in the case of the single maximum fabric is shifted downstream, but maintained in the simulation with single fabric in 2-D, while it disappears in 3-D.This leads to a different structure of age distribution in our main region of interest between the points V and H. Also, the horizontal velocity for the single vertical maximum scenario produces a completely different qualitative distribution from the 3-D case.Higher velocities are observed much further inside the valley (larger distances along the flow line) with the single fabric, while the simulations with the isotropic fabric are fairly similar (at least for distances beyond the outer first kilometer).

Discussion
We presented simulations with an isotropic and two different ad hoc prescribed fabric distributions: a depth-dependent distribution (evolving from isotropic at the surface to a vertical single maximum at one-third of the flow depth) and a uniform single vertical maximum fabric.The latter leads to a qualitatively correct reproduction of observed velocities at the BIA.
We get a similar picture from a comparison of the computed steady-state surface mass balance from Eq. ( 11) with the measured surface mass balance, shown in Table 3.Even   thickness in our simulations is neglected by the fixed surface, Table 3 makes it clear that the steady state solutions computed with isotropic as well as depth-dependent fabric produce surface mass balance estimates exceeding measured surface mass balance values at some places by more than an order of magnitude.Despite not obtaining a perfect match with observed data, in general the prescribed single vertical maximum fabric is the only tested fabric that resolves the paradox of not only lowering the average speed but also more strongly reducing the average vertical component of the velocity than the horizontal (see Fig. 3).Although the complexity of the flow regime prohibits an interpretation of the anisotropic model in detail, in a general way this can be explained because the ice becomes stiffer for vertical stretching or compression than shear.Since the anisotropic nature of ice introduces very strong differences in ice viscosity relative to the c axis of the ice crystal, introducing anisotropic fabric provides a plausible solution to the need to reduce flow speed in horizontal and vertical direction in a nonisotropic way.Other possible factors that would locally alter the ice viscosity, such as water content, ice impurities or damage, tend to soften the ice (Paterson, 1994) and increase the velocity, and consequently cannot explain the deviation from modeled isotropic velocity with respect to the observed at SBB. Errors introduced by the digital elevation models of either the bedrock or the free surface certainly have an impact on the flow velocity patterns.Nevertheless, noise in surface shape, by the incompressibility constraint of the fluid, would manifest in patchy flow patterns (see, e.g., Zwinger and Moore, 2009) rather than in the observed systematic deviation in direction and flow speeds.None of the 3-D simulations showed such behavior, leading to the conclusion that errors in the measured surface are not responsible for the mismatch of the velocities for the cases of isotropic and depth-dependent anisotropic fabrics.Since The Cryosphere, 8, 607-621, 2014 www.the-cryosphere.net/8/607/2014/ the blue ice areas are unusual in that the annual layers are tilted (varying from sub-horizontal near the equilibrium line to near vertical at the closed end of the ablation area), the fabric evolves over time and space.The full evolution of the fabric would have to take into account the non-steady-state evolution of the SBB valley since the LGM.During the LGM most of the valley was an accumulation region (Hättestrand and Johansen, 2005), and the blue ice was only formed as the surrounding ice sheet surface lowered and the mountain geometry created conditions suitable for effective removal of surface snow and firn (Malm, 2012).We could not simulate such a dynamic evolution of the valley given the considerable uncertainty of the date and timing of the ice sheet change that then drives the blue ice area evolution, and of course such a long time prognostic run would be prohibitively expensive in computer resources.The flow-line case suggested that both a depth varying fabric and a pure single maximum fabric would result in improvements to fit to observed velocities.In the three-dimensional simulation the depth-dependent fabric prescription (case (ii)) does not provide a satisfactory fit, while the single vertical maximum prescription (case (iii)) performs much better.The single vertical maximum fabric assumption is certainly too extreme a scenario and consequently leads to too low surface velocities.The real fabric distribution will be spatially varying, and our cases (i)-(iii) just provide a spectrum of the most extreme configurations.Nevertheless, the clear trend towards reducing the value and adjusting the directions of the horizontal velocity components while still matching the vertical ones (reflected in matching SMB numbers) give reason to claim that a pronounced fabric at the BIA of Scharffenbergbotnen remains the most likely explanation to correct the significant discrepancy between with isotropic rheologies computed and measured results.
From the comparison of the 2-D with the 3-D results for the age and the velocity in Fig. 10, it is clear that due to the complexity of the flow field, a full 3-D simulation is needed to capture all features.This is confirmed by the twisted distribution of the streamlines depicted in Fig. 9.
Our modeling results, in particular the difference between the prescribed depth-dependent and single vertical maximum fabric, suggest the presence of a strongly anisotropic ice fabric in SBB not only towards the bottom of the ice column but also near the surface.Considering that the ice fabric at a particular position depends not only on the local strain rates but also on the history of previous strain rates, our results indicate that the strongly anisotropic ice that typically appears in subsurface ice has been advected towards the surface of the BIA.
Strong fabric development has been observed in ice cores close to the surface (Svensson et al., 2007).Radar studies also show that the presence of strong fabric development near the surface is widespread in Antarctica (Matsuoka et al., 2012).The existence of such a strong anisotropy in ice fabric is explained by previous modeling results, which show that under compression, extension and shear, flowinduced anisotropy is developed in a fraction of the advection timescale.Such timescales are in SBB of the order or tens of thousands of years, defined by the ice thickness divided by vertical velocity (Martín and Gudmundsson, 2012).We have not modeled the evolution of the ice fabric here, mainly because of the lack of data, but our results support the idea that fabric evolves as rapidly in a BIA as in other areas with different ice-flow characteristics.This could be tested by collection of samples in the blue ice field; however this has never been done at SBB. Indeed such measurements should be made if flow modeling is to be used to help interpret horizontal ice cores.

Conclusions
We investigated the ice dynamics of SBB using a full-Stokes flow model as a tool to aid the interpretation of ice collected for paleoclimate research.By the time independence of the Stokes Eqs. ( 1) and ( 2), all components of the present-day velocity field are determined by the present-day geometry and current viscosity distribution within the ice, independently of the thinning history of the glacier.However, we cannot provide a transient forcing of the free surface and the fabric evolution since the LGM, which is a shortcoming when comparing model results to sampled ages and measured surface mass balances, as these properties are directly influenced by transient surface thinning.
We find that using standard isotropic ice fabric results in too high horizontal surface velocities and -in relation to the horizontal velocity scales -in too low vertical velocities.This makes for a very poor match with observations of surface age and motion.While it is possible that a temporal nonlinear thinning of the surface inside the valley since the LGM may contribute to this discrepancy, we find deviations by orders of magnitude between measured and computed surface mass balances for isotropic ice.Attempts to vary viscosity using plausible ranges of geothermal heat flux and surface temperature changes since the LGM could not resolve this issue.Within the assumption of steady-state dynamics, all other reasons for a drastic change of the velocity field with respect to the isotropic rheology solution may be dismissed as implausible; hence the only remaining explanation for the observed velocities is the existence of pronounced anisotropic fabric.
Confining our approach to look at extreme configurations of the rheology for a fixed geometry in SBB, we obtained the best match with the observed velocity field with a fabric showing single vertical maximum.As anisotropy appears to be a key factor in explaining the observed flow and agedepth distribution at SBB (and probably also in other blue ice areas), we suggest to investigate future ice samples not only with respect to chemistry but also with respect to the orientation of fabric.
www.the-cryosphere.net/8/607/2014/The Cryosphere, 8, 607-621, 2014 The remaining discrepancies concerning the match between the few reliable observed and computed ages at the surface certainly can be linked to the transient behavior of the age-depth profile.Even with a fixed geometry and prescribed fabric and temperature distribution (which is equivalent to a given viscosity), no steady state of the age profile can be obtained.Our simulated surface age in the inner part of the valley in any configuration is close to the integration time, and the maximum observed age at the surface in the innermost blue ice area is about 15 ka.This suggests that the BIA was formed around that time, which is plausible given the surface lowering that would have occurred in the ice sheet during the deglaciation.This is the first time that the importance of anisotropy in the ice dynamics of a blue ice area has been assessed.We show that the understanding of ice anisotropy is key in order to decipher the paleoclimatic records of blue ice areas.

Fig. 1 .
Fig. 1.Satellite pictures indicating (a) the position of Dronning Maud Land (credits Wikimedia Commons), with the SBB area indicated by the red dot, and (b) of the Scharffenbergbotnen (SBB) valley (NASA public imagery) and (c) a sketch of the domain augmented with observed stake and GPS velocities shown as arrows annotated by the stake reference code (red are GPS measurements, blue are theodolites) with corresponding stake positions as star and the color according to the velocity legend.GPS velocities are based on surveys using Svea station (large yellow circle) as the reference point.Geological and geographical features: brown rocks, dark-gray moraines, light-gray blue ice, orange lines distinctive features within the blue ice, long ice core in blue, position of the vertical transect for the flow-line model (red) within the area of SBB.All coordinates in this figure and also the following figures are defined by WGS84 UTM zone 29C.

Fig. 2 .
Fig.2.Sketch of the Scharffenbergbotnen (SBB) valley augmented with measured surface ages.Geological and geographical features are the same as in Fig.1.The inset on the upper left side is a close-up of the area around the horizontal ice core.Age measurements are shown as pentagons with the coloring according to the age legend with a distinction between reliable data (filled pentagons with black outlines) and data points with larger errors (empty pentagons with outlines in the color corresponding to their value) -see Table1.

Fig. 3 .
Fig.3.Summary statistics showing how component magnitudes of the velocity vector u = (u, v, w) T averaged over the inner BIA field are affected by changes in simulation parameters relative to the values in a standard simulation.Values are relative to simulations using an isotropic flow law (5); enhancement factor, E = 1; surface temperature, T zs = −20 • C; and geothermal heat flux, q geo = 60 mW m −2 (leftmost four columns).The other column sets represent results from simulations with a single parameter varying from this standard configuration as denoted by labels beneath the columns.BC represents a simulation with no inflow, u = 0, on the otherwise free inflow boundary of the domain.The rightmost set contains the results for the single maximum fabric simulation, and the inset emphasizes differences in horizontal and vertical components between this case and the isotropic case with the reduced enhancement factor (E = 0.1).

Fig. 4 .
Fig. 4. Horizontal velocity distribution (scale bar) obtained with the isotropic Glen's flow law compared with measured velocities (thick arrows; see Table2for values).Discrepancies in both direction and magnitude (with observed velocities being almost an order of magnitude smaller) can be easily seen.

Fig. 5 .
Fig. 5. Age (dots -scale bar) distribution obtained with the isotropic Glen's flow law and a 15 kyr integration of the age equation compared to measured (pentagons).Solid lines enclose BIAs.

Fig. 6 .Fig. 6 .Fig. 7 .
Fig. 6.Comparison of age (left column) and velocity (right column) values obtained on the flow-line (red line shown in Fig. 1) with the three-dimensional model using (a) isotropic Glen's flow law, (b) depth-dependent fabric and (c) single maximum fabric distribution (mind: one order of magnitude lower velocities).Age values are computed applying an integration over 15 kyr.The two vertical lines mark the intersection with the horizontal ice core at points V and H. 31 Fig. 6.Comparison of age (left column) and velocity (right column) values obtained on the flow line (red line shown in Fig. 1) with the three-dimensional model using (a) isotropic Glen's flow law, (b) depth-dependent fabric and (c) single maximum fabric distribution (note: one order of magnitude lower velocities).Age values are computed applying an integration over 15 kyr.The two vertical lines mark the intersection with the horizontal ice core at points V and H.

Fig. 7 .Fig. 8 .Fig. 8 .
Fig. 7. Surface velocity (vectors) for the anisotropic flow law using (a) a depth-dependent fabric and (b) a single maximum fabric distribution.The surface velocities in the depth-dependent anisotropic case (a) are not significantly altered with respect to the isotropic case.However, applying a single maximum fabric reduces velocities (even reducing them too much) and produces better agreement in direction.

Fig. 9 .
Fig. 9. Streamlines obtained with the backward Runge-Kutta method implemented in ParaView starting from the area where the ice cores have been sampled, colored by integration time given in years for (a) isotropic and (b) single maximum fabric distribution.The depth-dependent anisotropic fabric (not shown) is almost the same as the isotropic case and is characterized by stream lines originating in the area between the inner and the outer BIA and at the southern flanks.In contrast, the single maximum fabric shifts the origin of the stream lines towards the outer parts of the valley and significantly increases the travel time of the fluid particles.

Fig. 10 .Fig. 10 .
Fig. 10.Ages (left column) and absolute horizontal velocity (right column) of the 2-D flow-line model using (a) isotropic Glen's flow law, (b) single maximum fabric distribution (mind the different velocity scale).The two vertical lines mark the intersection with the horizontal ice core at points V and H.

Table 1 .
Measured calibrated ages at the surface including error in comparison to computed ages for the three different prescribed isotropic (iso), depth-dependent (depth) and the single vertical maximum (single) fabrics.The most reliable measurement points, V and H, are indicated in bold.Coordinate positions are given in UTM 29C.

Table 2 .
Measured surface velocities.Point sites within the inner BIA are indicated with bold letters.Coordinates are in UTM 29 C.

Table 3 .
Measured (Sinisalo et al., 2003)and computed steady-state surface mass balance in mm a −1 w.e. and comparison between measured and computed surface velocities for the isotropic (iso), depth-dependent (depth) and the single vertical maximum (single) fabric.Compare stake positions with Fig.1.