A full Stokes ice flow model for the vicinity of Dome Fuji, Antarctica, with induced anisotropy and fabric evolution

A three-dimensional, thermo-mechanically coupled ice flow model with induced anisotropy has been applied to a ∼200× 200 km domain around the Dome Fuji drill site, Antarctica. The model (“Elmer/Ice”) is based on the open-source multi-physics package Elmer ( http://www. csc.fi/elmer/ ) and solves the full Stokes equations. Flowinduced anisotropy in ice is accounted for by an implementation of theContinuum-mechanical, AnisotropicFlow model, based on an anisotropic Flow Enhancement factor (“CAFFE model”). Steady-state simulations for present-day climate conditions are conducted. The main findings are: (i) the flow regime at Dome Fuji is a complex superposition of vertical compression, horizontal extension and bed-parallel shear; (ii) for an assumed geothermal heat flux of 60 mW m −2 the basal temperature at Dome Fuji reaches the pressure melting point and the basal melting rate is ∼0.35 mm a−1; (iii) in agreement with observational data, the fabric shows a strong single maximum at Dome Fuji, and the age of the ice is decreased compared to an isotropic scenario; (iv) as a consequence of spatially variable basal melting conditions, the basal age tends to be smaller where the ice is thicker and larger where the ice is thinner. The latter result is of great relevance for the consideration of a future drill site in the area. Correspondence to: H. Seddik (hakime@lowtem.hokudai.ac.jp)


Introduction
Dome Fuji, situated at 77 • 19 01 S, 39 • 42 12 E, is the second-highest summit of the Antarctic ice sheet with an altitude of 3810 m a.m.s.l. (above mean sea level). Owing to its location on the Antarctic plateau (see Fig. 1) and the high elevation, it is one of the coldest places on Earth. Temperatures rarely rise above −30 • C in summer and can drop to −80 • C in winter. The annual average air temperature is −54.3 • C (Watanabe et al., 2003). The climate is that of a cold desert, with very dry conditions and an annual precipitation of about 27 mm water equiv., which falls entirely as snow (Kameda et al., 2008).
Deep ice core drilling at Dome Fuji was started in August 1995, and in December 1996 a depth of 2503 m was reached. The core was recently dated by an innovative method based on the O 2 /N 2 ratio of entrapped air (Kawamura et al., 2007), and covers a period back to ∼ 340 ka. A second deep core was started in 2003. Drilling was carried out during four subsequent austral summers from 2003/2004 until 2006/2007, and by then a depth of 3035.22 m was reached. The drill did not hit the bedrock, but rock particles and refrozen water have been found in the deepest ice, indicating that the bedrock is close to the bottom of the borehole. This core was dated preliminarily by synchronising the δ 18 O profile with that of the Dome C ice core data from the European Project for Ice Coring in Antarctica (EPICA), for which an age scale was already constructed (Parrenin et al., 2007a). As a result, the deepest parts of the second Dome Fuji ice core are estimated to be approximately 720 ka old, thus greatly extending the climatic record of the first core (Motoyama, 2007; Dome Fuji ice core project members, 2011).
Published by Copernicus Publications on behalf of the European Geosciences Union. The surface and bedrock topographies within 30 km from the Dome Fuji drill site were measured with high spatial resolution and resampled to a 0.5-km grid (Fujita et al., 1999;Watanabe et al., 2003;S. Fujita, personal communication, 2006). According to these data, the surface topography is smooth and flat, surface elevations 30 km away from the drill site being only 7 to 30 m lower than at the drill site. The actual summit is located ∼ 10 km west-north-west of the drill site and is ∼ 1 m higher. The bedrock topography exhibits a pronounced subglacial trench running mainly from the northeast to the south-west. Higher ground dominates to the southeast and north-west of the drill site, which itself is situated on a flank. The difference between highest and lowest bedrock elevations in the area is more than 1000 m.
The fabric (orientation distribution of the c-axes of the ice grains) of ice cores drilled near topographic domes typically evolves with depth from near-isotropic conditions close to the surface towards a single maximum. This was observed for the GRIP ice core in Greenland (Thorsteinsson et al., 1997), for the EPICA Dome C ice core in Antarctica (Wang et al., 2003) and also for the first Dome Fuji ice core, for which Azuma et al. (1999) report a strong vertical single maximum at 2500 m depth. Recently, A. Miyamoto (personal communication, 2009) extended the fabric measurements to 3024.5 m depth for the second Dome Fuji ice core. Results show that the strong single maximum remains essentially intact down to ∼ 2960 m. However, its direction is increasingly off the vertical, which coincides with layer inclinations up to 50 • (Dome Fuji ice core project members, 2011). Below ∼ 2960 m depth, the grain size increases significantly (probably due to migration recrystallisation in the warm, near-basal part of the core), and thin-sections of the ice core contain too few grains to allow a reasonable determination of the fabric.
Many models have been proposed to describe the anisotropy of polar ice. On the one end of the range in complexity, a simple flow enhancement factor is introduced in an ad-hoc fashion as a multiplier of the isotropic ice fluidity in order to account for anisotropy and/or impurities. This is done in most current ice sheet models, often without explicitly mentioning anisotropy (e.g. Saito and Abe-Ouchi, 2004;Greve, 2005;Huybrechts et al., 2007). In macroscopic, phenomenological models, an anisotropic macroscopic formulation for the flow law of the polycrystal is postulated. To be usable, the rheological parameters that enter this law must be evaluated as functions of the anisotropic fabric (e.g. Morland andStaroszczyk, 1998, 2003;Gillet-Chaulet et al., 2005. The concept of homogenisation models, also called micro-macro models, is to derive the polycrystalline behaviour at the level of individual crystals and the fabric (e.g. Lliboutry, 1993;Azuma, 1994;Castelnau et al., 1996Castelnau et al., , 1998Svendsen and Hutter, 1996;Gödert and Hutter, 1998;Ktitarev et al., 2002;Thorsteinsson, 2002). As for the "highend" complexity, full-field models solve the Stokes equation for ice flow properly by decomposing the polycrystal into many elements, which makes it possible to infer the stress and strain-rate heterogeneities at the microscopic scale (e.g. Meyssonnier and Philip, 2000;Mansuy et al., 2002;Lebensohn et al., 2004a,b). A very comprehensive, up-todate overview of these different types of models is given by Gagliardini et al. (2009). However, the more sophisticated models are usually too complex and computationally timeconsuming to be included readily in a model of macroscopic ice flow.
In this study, the Continuum-mechanical, Anisotropic Flow model, based on an anisotropic Flow Enhancement factor ("CAFFE model"; Placidi et al., 2010), is used. The CAFFE model belongs to the class of macroscopic models. It is a good compromise between physical adequateness and simplicity, and thus it is well suited to be embedded in flow models of ice sheets and glaciers. Since in the vicinity of topographic domes of ice sheets the widely used shallow ice approximation (Hutter, 1983;Morland, 1984) is not valid, we use the full Stokes model Elmer/Ice (e.g., Zwinger et al., 2007) instead, which is applied to an extended ∼ 200×200 km domain around Dome Fuji. The domain and the model equations, including the treatment of flow-induced anisotropy by the CAFFE model, are explained in Sect. 2. The numerical model itself is briefly introduced in Sect. 3. Section 4 constitutes the main focus of the study, where the simulation results for the ice flow, the anisotropic fabric, the temperature and the age are presented and discussed. Section 5 provides a discussion of some shortcomings of the chosen model approach, comparisons with available data and consequences for the consideration of another drill site in the region.

Coordinate system and domain
The vicinity of Dome Fuji is projected to a polar stereographic map with standard parallel at 77 • 19 01 S (the latitude of Dome Fuji) and central meridian at 0 • E. In the stereographic plane, we set up a Cartesian coordinate system with the south pole as origin. The horizontal plane is spanned by x and y, aligned with the natural x-and y-directions in Fig. 1, and z points upward. The associated unit vectors are e x , e y and e z , respectively. We define a quadratic domain of the size 203.092 × 203.092 km around Dome Fuji ( Fig. 1), which corresponds to a 200 × 200 km domain in the stereographic plane with standard parallel at 71 • S used by the Antarctica version of the large-scale ice sheet model SICOPOLIS (Iizuka et al., 2010). Within 30 km from Dome Fuji, the surface and bedrock topographies are given by the high-resolution data mentioned in the introduction. Outside this area, they are defined by the 20-km SICOPOLIS grid, based on resamplings of the RAM-PDEM V2 (Liu et al., 2001) and BEDMAP (Lythe et al., 2000) data-sets, respectively. We do not consider any evolution of the surface and bedrock topographies, but only deal with the present-day geometry. The topographies are shown in Figs. 2 and 3.

Field equations
Since ice is an (almost) incompressible material, conservation of mass entails that the velocity field (vector v) is solenoidal, div v = 0. (1) Further, the acceleration (inertia force) is negligible, so that the equation of motion is given by the incompressible Stokes equation (2) (e.g., Greve and Blatter, 2009), where p is the pressure, η the viscosity, ρ the ice density and g = −ge z the gravitational acceleration vector pointing downward. The viscosity is described by the CAFFE flow law, where d = 1 2 tr (D 2 ) is the effective strain rate, D = sym L = 1 2 (L + L T ) the strain-rate tensor (symmetric part of the velocity gradient L = grad v), n the power-law exponent, T = T −T m the temperature relative to pressure melting (T : absolute temperature, T m = T 0 − βp: pressure melting point, T 0 : melting point at low pressure, β: Clausius-Clapeyron constant), A(T ) the rate factor, A the polycrystal deformability andÊ(A) the anisotropic enhancement factor (see below, Sect. 2.3). The rate factor is expressed by the Arrhenius law where A 0 is the pre-exponential constant, Q the activation energy and R the universal gas constant. All parameters are given in Table 1. The temperature equation follows from the general balance equation of internal energy and reads Orientation transition parameter, ι 0.6 Length of year, 1 a 31 556 926 s (e.g., Greve and Blatter, 2009), where κ and c are the heat conductivity and specific heat of ice, respectively (Table 1). The age equation states that the age of ice, A, is equal to the elapsed time since an ice particle settled on the surface as a snowflake. This yields in a spatially fixed frame of reference the expression which describes a purely advective transport.

Boundary conditions
The ice surface, described by the function z = h(x,y), is assumed to be stress-free (atmospheric pressure and wind stress neglected), where t is the Cauchy stress tensor and n the outer unit normal vector. The surface temperature is prescribed according to ( Table 1), and the age of ice at the surface is set to zero, At the base, described by the function z = b(x,y), no-slip conditions are assumed, If the base is cold (temperature below the pressure melting point), the energy jump condition yields a Neumann-type boundary condition for the basal temperature, (e.g., Greve and Blatter, 2009), where the geothermal heat flux q ⊥ geo is prescribed (Table 1). If the base is temperate (temperature at the pressure melting point), the energy jump condition determines the basal melting rate, where L is the latent heat of ice. At the lateral boundaries (L.B.) of the domain, which are ∼ 100 km (∼ 30 times the thickness) away from the summit position, shallow ice stresses are prescribed, and the horizontal temperature gradient is assumed to vanish, (zero-flux condition), in order to minimise the influence on the temperature field in the domain.

CAFFE model
In order to take into account flow-induced anisotropy, we use the Continuum-mechanical, Anisotropic Flow model, based on an anisotropic Flow Enhancement factor ("CAFFE model"), which is explained in detail in the study by Placidi et al. (2010), and was already applied by Seddik et al. (2008),  and Bargmann et al. (2011) to the ice column of the EDML ice core in the Atlantic sector of East Antarctica. The polycrystal deformability A is defined as where S 2 is the unit sphere, S is the deviatoric stress tensor, n is the orientation (normal unit vector of the basal plane, direction of the c-axis) and f * (n) is the orientation distribution function (ODF). The latter is defined as where ρ * (n) is the orientation mass density (OMD), that is, the mass per volume and orientation. By introducing the second-and fourth-order orientation tensors, the polycrystal deformability (Eq. 15) can be expressed as The anisotropic enhancement factorÊ(A) depends on the polycrystal deformability by the smooth, monotonically increasing function with the two parameters E max = 10 (maximum softening) and E min = 0.1 (maximum hardening). The function is illustrated in Fig. 4. The evolution equation for the OMD is the orientation mass balance In this equation, div is the normal, three-dimensional divergence operator, div S 2 the divergence operator on the unit sphere, u * (n) the orientation transition rate, q * (n) the orientation flux and * (n) the orientation production rate. Physically, the orientation transition rate corresponds to grain rotation, the orientation flux to rotation recrystallisation (polygonisation), and the orientation production rate to migration recrystallisation. Like in the study by Seddik et al. (2008), we do not take into account recrystallisation processes here, and thus set q * (n) = 0, * (n) = 0. This leaves the orientation transition rate to be determined, for which the constitutive equation is employed. The parameter ι is a positive constant. The additional term W · n with the spin tensor W = skw L = 1 2 (L − L T ) (skew-symmetric part of the velocity gradient L) describes the contribution of local rigid-body rotations. For a spatially three-dimensional problem, the evolution equation (Eq. 20) is six-dimensional (three spatial directions, time, two orientation components), which renders its numerical solution on a reasonably fine grid unfeasible. Therefore, we resolve to a simplified formulation with the orientation tensors introduced above. By differentiating Eq. (17) 1 with respect to time and inserting the constitutive equation (Eq. 21), one obtains an evolution equation for the secondorder orientation tensor, (for the detailed derivation see Seddik, 2008). This equation contains the unknown fourth-order orientation tensor a (4) on the right-hand side. We formulate a closure condition of the form that is, we relate a (4) by some non-linear function F IBOF to a (2) . Specifically, following Gillet-Chaulet (2006) and , we adopt the invariant-based optimal fitting closure approximation (IBOF) proposed by Chung and Kwon (2002) because of its physical accuracy and numerical stability compared to previously proposed closure functions. The exact form of the very lengthy function F IBOF shall not be given here; we refer the reader to ,  or Seddik (2008) instead. As for the orientation transition parameter ι, Placidi (2004) showed that the fabrics in the upper 2000 m of the GRIP ice core can be best explained by the value ι ≈ 0.4, while the study on the East Antarctic EDML ice core by Seddik et al. (2008) provided a best fit between modelled and measured fabrics for ι = 0.6. Since the region of the latter study is close to our one (both are located in Dronning Maud Land, East Antarctica) and the conditions are similar, we follow Seddik et al. (2008) and adopt the value ι = 0.6. The boundary condition for Eq. (22) is isotropy at the ice surface, that is, where I is the unit tensor.

Finite element model Elmer/Ice
The Using the commercial pre-processing tool Gambit, a mesh of the computational domain defined in Sect. 2.1 has been created. The mesh has a rather coarse resolution of 7 km near the boundaries, which refines to 0.5 km towards the centre of the domain, located very close (∼ 0.23 km) to the position of Dome Fuji. An inner square domain (6 × 6 km) is meshed with mapped elements with a resolution of 0.5 km. This domain is surrounded by a second square domain (60 × 60 km) meshed with mapped elements using a 6-km resolution. The connection between those two square sub-domains is performed with a refining resolution towards the inner square. In order to connect the outer square to the domain boundaries, four sub-domains have been created and meshed with paved elements. The vertical direction is discretised by 20 equidistant layers. The resulting mesh contains 5406 bilinear quadrilateral elements, 30 860 trilinear brick elements and a total of 33 642 nodes (Fig. 5).
The numerical solution technique is similar to the one explained in the study by Zwinger et al. (2007) in the context of an application of Elmer/Ice to a small crater glacier. In particular, we only consider steady-state scenarios with presentday climatic forcing and fixed geometry. Differences are that we are concerned with incompressible ice instead of compressible firn, and that we take into account flow-induced anisotropy by the CAFFE model. In detail, we first compute a steady-state solution of the velocity and temperature fields for a linear rheology (power-law exponent n = 1 instead of 3) and isotropic conditions (enhancement factor uniformly set to E = 1). In the second step, we use the results as initial conditions for an isotropic steady-state run with n = 3. The velocity and temperature fields produced by this run serve as initial conditions for the final, fully coupled simulation, which is integrated over 1 Ma with temporally constant boundary conditions in order to obtain steady-state solutions of the velocity, temperature, age and anisotropic fabric. The model employs a time step of 100 a during the first 100 ka and a time step of 200 a during the subsequent 900 ka. As for the solution of the advective fabric evolution equation (Eq. 22), we add some numerical diffusion N of the form N = ν hor ∂ 2 a (2) ∂x 2 + ∂ 2 a (2) ∂y 2 + ν vert ∂ 2 a (2) ∂z 2 (25) on the right hand side in order to ensure stable solutions. The quantities ν hor and ν vert are the horizontal and vertical diffusivities, respectively. Additional boundary conditions at the ice base are now required, and we choose the Neumann-type conditions The numbers 1/3 in Eq. (26) 1,2 result from the surface boundary condition (Eq. 24), the numbers 0.9 and 0.05 are rough a priori guesses of the basal values of a zz and a xx/yy at Dome Fuji, respectively, and H DF is the ice thickness at Dome Fuji. Thus, Eq. (26) 1,2 represent estimates of the basal gradient of a basal values and the assumption of linear profiles. For the non-diagonal elements, a zero-gradient condition is assumed (Eq. 26 3 ).

Anisotropic steady-state simulation
We will now discuss the results of a present-day steadystate run for the 203.092×203.092 km domain around Dome Fuji as described in Sects. 2, 3 and Table 1. Values of the geothermal heat flux exceeding q ⊥ geo ∼ > 50 mW m −2 result in a temperate base below the Dome Fuji drill site. We choose q ⊥ geo = 60 mW m −2 (see Table 1), because the additional heat flux produces a basal melting rate of ∼ 0.35 mm a −1 at the Dome Fuji drill site, which leads to a basal age of ∼ 635 ka (Fig. 15). This is close enough to the estimated basal age of 720 ka for the real Dome Fuji core reported by Motoyama (2007) and the Dome Fuji ice core project members (2011). In order to prevent the basal age from becoming infinite in regions where the basal temperature is below the pressure melting point and the ice is frozen to the ground, a basalmelting offset of 0.1 mm a −1 is applied there. The numerical diffusivities in Eq. (25) are chosen as ν hor = 10 −5 m 2 s −1 and ν vert = 10 −6 m 2 s −1 , which will be justified below ( Fig. 9 and associated main text). Figure 6 shows the simulated surface velocity for the entire domain. The ice flow points away from the summit, follows largely the direction of the steepest surface slope despite the pronounced relief of the bedrock (Figs. 2, 3) and reaches maximum speeds of ∼ 1.6 m a −1 around the margin of the  . The very small strain rates at the Dome Fuji drill site shown in Fig. 7 (of the order of 10 −5 a −1 ) result from the fact that the drill site is only ∼ 10 km away from the actual summit position. As expected, the flow regime at the drill site is a complex superposition of vertical compression (D zz ), horizontal extension (mainly D xx ) and bed-parallel shear (D xz , D yz ) (e.g., Raymond, 1983). Vertical compression and horizontal extension dominate in the upper parts of the ice column, whereas bed-parallel shear becomes more important in the lower parts.
The largest eigenvalue a 3 (that corresponds approximately to the vertical direction) of the second-order orientation tensor a (2) near the ice base (at 95 % depth), which is shown in Fig. 8 found in the immediate vicinity of the topographic summit (∼ 10 km west-north-west of the drill site; see the introduction), and indicate the presence of a strong single maximum fabric. This is a consequence of the regime of vertical compression which dominates there. Further away, in the areas where bed-parallel shear prevails, lower values between ∼ 0.55 and 0.65 are found, which correspond to weakly to moderately developed fabrics. Directly at the drill site, Fig. 9 shows the anisotropic fabric expressed as the smallest (a 1 ), middle (a 2 ) and largest (a 3 ) eigenvalue of a (2) over depth. The anisotropy manifests itself by values of a 3 larger than and values of a 1 and a 2 smaller than 1/3. A strong single maximum fabric is found at greater depths. Since a 1 < a 2 , the maximum is slightly elongated, which is consistent with the the dominance of the horizontal extension component D xx over D yy (Fig. 7).
As for the values of the numerical diffusion, the factor 10 difference between ν hor (standard value 10 −5 m 2 s −1 ) and ν vert (standard value 10 −6 m 2 s −1 ) reflects the shallowness of the model domain and ensures that the numerical stabilisation is efficient in both the horizontal and the vertical direction. Figure 9 also shows the profiles of a 3 for 10 times larger and smaller diffusivities, respectively. For 10 times larger diffusivities, the profile is essentially a straight line with a slope equal to the prescribed gradient at the base, which indicates that the influence of the numerical diffusion on the result is too large. For 10 times smaller diffusivities, the profile is very similar to that computed with the standard values down to ∼ 1500 m depth, but then shows a divergence and, particularly, a spurious, unrealistic decline further down. For this particular test, it was not possible to reach 1 Ma of sim- ulation time; the simulation becomes numerically unstable and stops at 330 ka. This observation emphasizes the need to use sufficiently high numerical diffusivities. The results obtained with the boundary conditions from Eq. (26) are also compared to results obtained with a zerogradient boundary condition for all elements of a (2) (blue line in Fig. 9). The differences of the values of a 3 are significant; while limited to ∼ 10% in the upper half of the core, they grow to ∼ 30% towards the base and lead to a weaker fabric. Hence, a reasonable a-priori guess of the basal gradient of a (2) (like in Eq. 26) is necessary in order to get a reasonable result. This situation, while unavoidable at the moment, is not entirely satisfactory and deserves further thoughts.
In an additional test, we have replaced the uniform ice thickness H DF in Eq. (26) by the local ice thickness H , which reduces the prescribed gradient of a (2) where the ice is thicker than at Dome Fuji and vice versa. In that case, the differences compared to the results shown in Figs. 8 and 9 are minimal, so that the different choices of the ice thickness do not matter. The simulated basal temperature relative to the pressure melting point is depicted in Fig. 10. Approximately 13 % of the basal ice in the domain is at the pressure melting point, including the Dome Fuji drill site itself. These areas are characterised by large ice thickness, whereas lower basal temperatures (down to ∼−10 • C) occur in areas where the ice is relatively thin due to higher bedrock elevations (Fig. 3). The temperature profile of the Dome Fuji drill site (Fig. 11) shows a continuous increase with depth and is only moderately concave, which results from the small contribution of advection to the total heat transport. The simulated profile agrees well with data measured in the field season 2006/2007 (only below 200 m depth; H. Motoyama, personal communication, 2010). It must be noted, however, that due to technical problems (no compensation circuit was available for the measurements with a Pt sensor) the data may suffer from a systematic error that cannot be assessed quantitatively.
The simulated vertical velocity profile is depicted in Fig. 12. The values vary from ∼ 41.1 mm a −1 at the surface to ∼ 0.35 mm a −1 at the base, which corresponds to the basal melting rate. The simulated profile is compared to a one-dimensional flow model devised by Parrenin et al. (2007b), where the surface accumulation, basal melting and velocity profile were reconstructed using age markers as constraints. The two vertical profiles are essentially identical in shape; however, the one-dimensional model shows significantly smaller values, 29.1mm a −1 at the surface. The  Motoyama, personal communication, 2010). latter is very close to the present-day surface accumulation rate of 27 mm water equiv. a −1 = 29.7 mm ice equiv. a −1 (Kameda et al., 2008), while the value obtained by our simulation is 38 % larger. This discrepancy is likely due to the fixed-topography, steady-state approach used in our simulation. The basal melting rate reconstructed by Parrenin et al. (2007b) is 0.011 mm a −1 , which is much smaller than our value of ∼ 0.35 mm a −1 . However, the value by Parrenin et al. (2007b) is very poorly constrained because their model is based on the first Dome Fuji ice core only, of which the bottom is 500 m above the bedrock. The authors assess the range of uncertainty of the basal melting rate as less than 0.4 mm a −1 , with which our value agrees. Figure 13 shows the simulated age at 95 % depth for the entire model domain. At the Dome Fuji drill site, its value is ∼ 350 ka and increases to ∼ 635 ka directly at the base. This is somewhat lower than the estimated real value of 720 ka (Motoyama, 2007). Of course, our results are biased by the simplifying steady-state assumption, and we do not aim at providing an accurate dating in terms of absolute values. However, the distribution of the age is less affected by this assumption and should be reasonable. The relation to the ice thickness (Fig. 3) is not evident at first glance. The scatter plots shown in Fig. 14  positive at 95 % depth and strongly negative at the base; that is, the basal age tends to decrease with increasing ice thickness (note that Fig. 14b shows only the results for warmbased areas because of the above-mentioned problem with ill-defined ages in case of a cold base). The reason for this behaviour is that basal melting increases with increasing ice thickness (∼ 0.35 mm a −1 at the Dome Fuji drill site), which melts away the oldest ice (e.g., Fahnestock et al., 2001).

Comparison with an isotropic simulation
In order to investigate the influence of the flow-induced anisotropic fabric, we have also conducted a present-day steady-state simulation with the fabric evolution "switched off" and the enhancement factor uniformly set to E = 1 instead. The resulting surface velocity distribution of this isotropic simulation turns out to be very similar to that of the anisotropic simulation shown in Fig. 6 slightly smaller in the isotropic case, in agreement with previous findings by Mangeney et al. (1996) and Pettit et al. (2007). This affects the age of the ice: the slightly smaller flow velocities lead to slower downward transport, so that the isotropic simulation produces generally larger ages than the anisotropic simulation. This is illustrated in Fig. 15, which depicts the age profiles at Dome Fuji for both simulations. At the base, a value of ∼ 660 ka from the isotropic simulation contrasts with ∼ 635 ka from the anisotropic simulation. It is also interesting to note that the slightly positive correlation between the age at 95 % depth and the ice thickness, as well as the strongly negative one between the basal age and the ice thickness, observed above (Fig. 14) are robust features and show up also in the isotropic simulation (Fig. 16). The slopes of the corresponding regression lines are somewhat larger in the isotropic case, which demonstrates the nonnegligible effect of anisotropy on age; however, basal melting limits efficiently the age of the basal ice for large thicknesses irrespective of the consideration or non-consideration of anisotropy.

Discussion
Some caveats to the results presented above must be noted. Since the Antarctic ice sheet has been subjected to largely varying climatic conditions in the past (glacial-interglacial cycles), our steady-state approach constitutes a major simplification. The topography in the vicinity of Dome Fuji depends crucially on the dynamical history of the Amery Ice Shelf and its drainage area. Saito and Abe-Ouchi (2010) show that, if the Amery Ice Shelf was entirely grounded during a prolonged period in the past, the topographical summit at Dome Fuji may have vanished completely, which would have exposed the area to normal simple-shear flow. Hence, it is possible that the dynamical regime in the vicinity of Dome Fuji experienced significant changes over time. A further source of uncertainty is the assumption of a spatially constant geothermal heat flux. Even for a rather small area like our ∼ 200 × 200 km domain this is not necessarily the case, and distributions of the basal temperature, melting rate and age would be affected by a considerable spatial variability.
In this context, it is interesting to note that the layer inclination of the Dome Fuji ice core shows sudden increases at several depths and reaches ∼ 50 • near the bottom (Dome Fuji ice core project members, 2011). This can be explained partly by the contrast between cold-based ice (no basal melting) and temperate-based ice (basal melting) in the immediate vicinity of the drill site found in our simulations (Fig. 10). However, the extreme layer inclination near the bottom and especially the discontinuous increase with depth probably require additional mechanisms. One possibility is that the spatial distribution of the geothermal heat flux has experienced abrupt changes at several occasions. Alternatively, folding (Thorsteinsson and Waddington, 2002) or episodic drainage of a nearby subglacial lake (Pattyn, 2008) may be responsible for the irregular layer inclination.
The simulated fabric shown in Fig. 9 agrees well with the data by Azuma et al. (1999) who report an essentially monotonic transition (with a few irregularities, though) from isotropic conditions at the surface to a strong vertical single maximum at 2500 m depth. As already mentioned in the introduction, the data by A. Miyamoto (personal communication, 2009) show an essentially intact, but increasingly offvertical single maximum at greater depths. The latter is not reproduced by our simulation and most likely related to the extreme layer inclinations discussed above. The observed behaviour of rapidly increasing grain size below ∼ 2960 m depth indicates the onset of migration recrystallisation due to the warm conditions in the near-basal part of the core. This is not accounted for in our study (see the discussion of the orientation mass balance; Eq. 20), and thus any effect on the fabric (which cannot be inferred from the data due to too large grains) is not reproduced.
All of these conjectured mechanisms can alter the distribution of the age of the ice in the region. Furthermore, the old ice predicted by our simulations for rather thin locations suffers from strong layer thinning, which is less pronounced for the younger ice at thicker locations. With regard to the idea of retrieving another deep ice core in the vicinity of Dome Fuji, this is important to consider because the temporal resolution of the oldest, near-basal ice would be quite poor. An area-wide radio-echo-sounding survey, possibly supplemented by melting probes at selected positions, would be desirable in order to obtain more information about the detailed conditions in the region.