Numerical simulations and observations of the role of katabatic winds in the creation and maintenance of Scharffenbergbotnen blue ice area, Antarctica

We model the role of katabatic winds in the formation and maintenance of a blue ice area in Scharffenbergbotnen valley, Antarctica, using the finite element code Elmer. The high-horizontal-resolution (50–200 m) numerical simulations of the local wind flow from katabatic wind fronts show high spatial variability in wind-impact patterns and good congruence between places with high near-surface wind speeds and the blue ice area. In addition we perform wind simulations on an altered glacier geometry that resembles the thicker ice cover at the Late Glacial Maximum (LGM). These simulations indicate that the pronounced spatial wind-impact patterns depend on present-day geometry and did not occur during the LGM. This leads to the conclusion that the formation of the inner blue ice area of the Scharffenbergbotnen valley started only after the lowering of the ice surface, i.e. after the LGM. Experiments with smoothed surface topography suggest that detailed positions of the high wind regions, and hence individual blue ice fields, may have varied as the ice sheet lowered. The simulation results obtained with the present-day geometry were fortuitously confirmed by the destruction of a field camp located in a high-wind-speed area and its subsequent redistribution to low-velocity areas. The experiments and the field observations are consistent with localized violent katabatic events rather than synoptic-scale storms, playing the dominant role in the formation and maintenance of this and perhaps many blue ice areas.


Introduction
A blue ice area (BIA) is a snow-free ablation zone where the surface ice reflects a blueish colour. These areas cover about 1 % of the entire surface area of Antarctica (Sinisalo and Moore, 2010). At relatively low-elevation BIAs, summer temperatures are sufficiently high to allow surface and internal melting, which is of great practical importance for water supply to field stations and landing strips for wheeled aircraft. The lower surface roughness of ice compared with snow-covered surfaces allows winds to preferentially remove any freshly fallen snow from them. At higher altitudes, typically above 1000 m a.s.l. (Bintanja, 1999), melt does not occur even in summertime, and BIAs are then maintained and created simply by the occurrence of strong winds over them. It is one of these high-elevation BIAs that is object of study in this paper. The snow-ice boundaries of BIAs appear to be remarkably stable on multi-year timescales (Giaever, 1969), which indicates that repeating meteorological conditions must play a role in maintaining them (Bintanja, 1999;Sinisalo et al., 2004;Sinisalo, 2007). Old ice can often be found at the surface of some BIAs. If a good understanding of the ice dynamics were available, then the easily accessible ancient ice would be of great interest for palaeoclimatic studies (Bintanja, 1999;Sinisalo et al., 2004;Moore et al., 2006;Sinisalo and Moore, 2010).
BIAs are hence typically located in areas where the local wind speed is occasionally high -either by association with Figure 1. A map of Antarctica is shown at the top right with the domain shown in main image in the boxed area. In the main image, the AWS6 and AWS9 locations are marked and the Scharffenbergbotnen area is enclosed by the box and shown in detail at the top left. Wind tails can be seen generally in a westerly direction from the nunataks, indicating winds from the east. The green triangle on the eastern border of the larger inner BIA shows where a field camp was located before it was destroyed by a storm on 2-3 January 2007, and the magenta triangle shows where one destroyed tent was found after the storm with debris scattered over a wide area of the snowcovered ground. Imagery courtesy of NASA. mountain turbulence or in regions of strong katabatic wind flow. Here we investigate to what extent the surrounding terrain contributes to a local increase in wind-impact velocity above the well-studied BIAs in the inner region of Scharffenbergbotnen (SBB) valley, Dronning Maud Land, Antarctica (located approximately at 74.56 • S, 11.05 • W), and compare the wind distribution patterns with the equilibrium surface mass balance distribution computed by Zwinger et al. (2014). SBB is located about 300 km inland, in the Heimefront mountain range of Dronning Maud Land (see Fig. 1). The dominant easterly wind direction can be clearly seen from the prominent wind drift tails behind exposed rocks. A geomorphological map of the SBB valley is presented in Fig. 2. The map shows three distinguishable BIAs, one on the north-western side at the entrance of the valley and two inside the SBB valley, separated by a moraine. The two BIAs inside the valley have been the subject of various studies of meteorological (Bintanja and Reijmer, 2001), glaciological (Sinisalo and Moore, 2010) and geomorphological nature (Hättestrand and Johansen, 2005).
Numerous previous studies of winds over Antarctica have adopted simplified equations and worked on large spatial domains, ignoring small-scale details and only resolving largescale flow patterns. A hydrostatic three-dimensional primitive equation model with a resolution of 20 km was presented by Bromwich et al. (1994) for simulations on winter katabatic winds crossing the Siple coast and the Ross ice shelf. Gallée and Schayes (1994) proposed a similar model with the highest resolution being 5 km. A relatively simple simula- Figure 2. The Scharffenbergbotnen area, with the location of the blue ice areas (purple), exposed rocks on top of hills and mountains (brown) and moraines (grey). Isolines of elevation are given in steps of 25 m. The inner BIAs are divided by a moraine into a northern and a southern BIA. Many of the other marked moraines are actually blue ice covered by rocks and much of the flatter area south and east of the inner blue ice areas is likely blue ice more or less covered by a thin layer of rocky material carried by ice flow from snow accumulation regions on the valley sides. Zwinger et al. (2014) show that ice flow to the inner northern BIA originates from snow accumulated on the ridge between it and the outer BIA, indicated by the long arrow in the figure. The green and magenta triangles represent the destroyed camp and location of debris (as in Fig. 1). Coordinates are in UTM 29 C system. tion of snow flow affected by katabatic winds was presented by Gallée (1998); with a horizontal resolution of 40 km, the model domain covered the entire Antarctic ice sheet. The study by Yu et al. (2005) on katabatic jumps, where an open channel katabatic flow rapidly transforms to a calm flow, was implemented on a uniform grid spacing of 3.5 km, 700 and 175 m along with an assumption of a two-dimensional katabatic flow and a rather smooth geometry. Skyllingstad (2003) in contrast used a large eddy simulation (LES) to investigate the small-scale dynamics of gravity-driven incompressible flows on a simple linearly decreasing smooth slope with a resolution of approximately 0.75 m. Axelsen and van Dop (2009) applied LES to simulations of steady-state quasi-onedimensional buoyancy-driven slope winds over glaciers.
We deploy variational multiscale (VMS) stabilization with the finite element (FE) method to study the local wind fields above the valley at SBB. To our knowledge this is the first high-resolution Navier-Stokes governed simulation of local wind fields above BIAs. The domain size used for our modelling of the complex terrain in SBB is 10 3 times larger than in Skyllingstad (2003), and our minimal horizontal resolution is h min ≈ 50 m.
First we present observational evidence that the mechanism responsible for creating and maintaining the inner BIA at SBB is to be found in easterly katabatic wind fronts in Sect. 2. This provides motivation for the following analy- sis. In Sect. 3 we present the governing equations including a brief introduction to the VMS turbulence model in Sect. "variational multiscale method". This will be followed in Sect. 4 by a discussion of the specific model set-up for the simulations of SBB area. Section 5 presents and discusses the results obtained for different flow profiles as well as altered geometries. The conclusions from these results are presented in Sect. 6.

Case study of observed violent wind event
In January 2007 a field camp located near the head of the Scharffenbergbotnen valley was destroyed by a strong storm (Fig. 3). There are two reasonably close automatic weather stations (AWSs), AWS6 and AWS9, in the region of Scharffenbergbotnen. AWS9 is located at the summit of the regional ice field at 75 • 00 09 S, 00 • 00 26 E and 2900 m a.s.l., while AWS6 is about 40 km north of the Scharffenbergbotnen valley at 74 • 28 53 S, 11 • 31 01 W and 1160 m a.s.l. (Bintanja, 1999;Van den Broeke et al., 2005); the location of the weather stations in relation to Scharffenbergbotnen is shown in Fig. 1. Figure 4 shows the wind direction and magnitude recorded by the two AWSs over the first 22 days of 2007.
The largest wind speeds occur on the 11th and 12th days and are consistent in direction at the two AWSs, with somewhat higher wind speeds occurring, as expected, at the lower AWS6 site. However, the wind speeds in the Scharffenbergbotnen valley during the 11th and 12th days were not very high; indeed field work progressed normally on both days, with no drifting snow, and the 12th day was notably sunny and pleasant. This was a synoptic-scale event that may be  classed as average conditions. However, the event of 2-3 January was extremely strong in Scharffenbergbotnen but not significant on a synoptic scale, with moderate winds and inconsistent directions at AWS6 and AWS9. Such relatively stable conditions on the higher ice sheet may actually favour the formation of a katabatic jet. Even 40 km north of the Scharffenbergbotnen valley the katabatic jet had dissipated sufficiently for there to be no trace of it in the AWS6 record. The debris from the 2-3 January storm was dispersed over much of the valley and provides traces for recirculation areas. The location where the field camp was originally placed and where a tent from the destroyed camp was found is shown in Fig. 1. This and many other items were all found on snow accumulation areas and not on the blue ice area. Existing snow and even hard firn was ripped from the BIA such that the extent of the BIA was increased by perhaps 30 % relative to that before the storm. This suggests that such storms are rather infrequent and that the BIA extent may be governed by rare stochastic events that remove the steadily accumulated snow and firn of several years. In what follows, we present a numerical model that studies the impact of a turbulent wall-bounded jet into the Scharffenbergbotnen valley that enables us to identify the areas of high wind speeds close to the ground within the valley where potentially snow can be removed.
The flow is assumed to be described by the incompressible Navier-Stokes equations where ρ stands for the constant density of the air, f for the gravitational body force, u denotes the velocity vector, p the isotropic pressure and µ is the dynamic viscosity . The numerical values of these parameters used in the simulation are presented in Table 1. The symmetric strain rate tensor, ε, is given by The set of Eqs. (1) and (2) is defined in the domain ⊂ R 3 for the time interval t ∈ (0, t end ).
We are well aware that katabatic winds are driven by the density differences between colder and warmer air and hence a detailed study would deserve a buoyant (hence compressible or at least Boussinesq-approximated) flow with its density coupled to the temperature field. Nevertheless, our study will show that over the last few hundred metres of the katabatic flow into the valley the inertia of the incoming jet is sufficient to produce wind-impact patterns that match the observed locations of the bare ice fields in the SBB valley well, and hence the assumption of incompressibility is appropriate for the BIA study. In Appendix A we present Eq. (1) in its non-dimensional form, accounting for changing density using a Boussinesq approximation. By an analysis of the non-dimensional groups we can demonstrate that in the case of a fully developed katabatic jet for the length scales that apply in our model, the inertial forces dominate the buoyancy forces. This shows the applicability of an incompressible model to this work.

Variational multiscale method
The modelling of turbulent flows is one of the most challenging areas in mathematics and computational fluid dynamics (CFD). Numerous ad hoc models to treat the turbulent nature of the resulting flow numerically have been proposed. The main idea of using these approaches is to dramatically reduce the number of nodes needed for solving turbulent CFD problems (Rautaheimo, 2001) in comparison to direct numerical solution. Sub-grid-scale turbulent models, such as LES models, resolve the large-scale turbulence and model the effects of sub-grid turbulence. In the case of LES the sub-grid model is obtained by spatial filtering, which imposes certain problems in particular with bounded flows.
We use the implicit Euler method for time discretization. In space, continuous piecewise isoparametric hexahedral elements are used for both the velocity and pressure.
The discretization in space is a weighted least squares stabilized Galerkin method, also known as the VMS method (Hughes, 1995;Johnson, 2006, 2007;Bazilevs et al., 2008). The VMS method is an alternative, yet similar, method to LES. The main difference between LES and VMS is that the latter projects the scales into the function spaces rather than filtering them by spatial averaging. Thus the VMS method avoids the problems of spatial filtering, spatial differentiation in bounded domains as well as with nonconstant filter widths and the choosing the correct boundary conditions for spatially averaging at large scales (John et al., 2010).
Given the initial value u 0 h , the solution (u n h , p n h ) for the time step t n = n t, n = 1, 2, . . ., N , is obtained from where (·, ·) denotes the L 2 ( ) inner product and (·, ·) E denotes the L 2 (E) inner product for the element E of the mesh. V 0 h and W h represent the finite element spaces. A constant time-step size of t = 0.1 s is applied in all simulations. The stabilization parameter is given by where h E represents the size of element E. Detailed literature on applications and the theory behind the VMS method can be found in Hughes (1995), Gravemeier (2003), Gravemeier et al. (2006) and Bazilevs et al. (2008). At each time step the nonlinear system (Eq. 3) is solved by a Picard iteration.
The linearized system is solved using the generalized minimal residual Krylov subspace method (Kelley, 1995) with an incomplete LU (lower upper) factorization of the system matrix deployed as a preconditioning step.

Simulation set-up for Scharffenbergbotnen valley
The VMS method used for turbulent computations is implemented within the open-source multi-physics FE-code Elmer (http://www.csc.fi/elmer). Due to a large domain size combined with a turbulent flow the problem has to be solved using high-performance computing on a parallel Linux cluster. Thus a sufficiently fine mesh in the areas of most interest can be achieved while avoiding excessive computations. As shown in Fig. 5, this footprint mesh is vertically extruded into 30 terrain-following layers up to an elevation of 3000 m. The first 10 mesh layers are concentrated into the lowest 10 % of the height, applying an exponential refinement of the layers towards the lower surface. Thereby a high vertical resolution (1-2 m) at the valley bottom to account for the turbulent boundary layer close to the surface is achieved. This results in a mesh with 373 620 nodes from which 358 382 (mainly hexahedral) linear bulk elements are constructed, which are bounded by 30 516 (mainly quadrilateral) linear boundary elements. This mesh is split into 64 partitions for parallel computation.
Our simulation set-ups resemble those from Malm (2012). However, while Malm (2012) deployed initial conditions utilizing a background wind field statistically obtained from ECWMF reanalysis data, we let the katabatic wind front impact into a region of stagnant air. No significant differences in wind-impact patterns between these background fields were found.
In addition to using present-day topography, two alternative meshes were constructed: a smoothed version of the initial DEM and one based on geomorphological evidence of the ice thickness at Late Glacial Maximum (LGM), i.e. a DEM with 200 m thicker ice inside and 100 m outside the valley (Näslund et al., 2000;Hättestrand and Johansen, 2005). To approximate thicker glacial ice conditions around the valley, all non-glaciated peaks outside the valley were also smoothed. The smoothing in both variants of the original DEM was achieved by replacing the elevation value of each point in the smoothed DEM by the average of all surrounding points contained within a circle of 250 m radius in the original DEM, similar to a simple moving average algorithm.
We also allow for variation in the vertical velocity profile of the katabatic wave front as it enters the eastern border of the domain. The reference simulation profile for the normal, inwards-pointing component of the velocity is denoted as nose shaped, which is defined with respect to the height above ground, H [m] as The alternative is a parabolic profile (the same as used by Malm, 2012): Both profiles reach from ground level to 100 m above ground. The third profile investigated is a stretched version of Eq. (5), noseL, which extends to 200 m above ground: All three profiles are shown as a function of H in Fig. 6. To avoid issues with backflow at the lateral boundaries, the inflow profile is damped to zero at a distance of 1250 m from the lateral (north and south) boundaries. Instabilities imposed by the pressure field at the start of the simulation were suppressed by allowing 100 s for the inflow profile to ramp up from resting air to peak values. On the southern and northern boundaries, which are aligned with the main flow direction, the normal component of the velocity is kept at 0. The same condition applies to the boundary with the upper atmosphere, where we constrain the vertical (z aligned) velocity component. The set of boundary conditions is completed at the valley surface, b(x, y), where we impose a no-slip condition, and a further free outflow condition imposed on the western confinement. In order to avoid instabilities, any backflow at this boundary is reset to 0, i.e.
u · n | west ≥ 0, where n is the outward-pointing normal of the boundary.

Results
With the set-up discussed in the previous section, we study the impact of wind on the valley. Our main focus is to study the BIA north of the moraine inside the valley, i.e. the larger one stretching from about 435 to 439 km easting and 1722.5 and 1724 km northing indicated as inner northern BIA in Fig. 2. A first approach is to test the original DEM with the noseshaped profile defined in Eq. (5) and compare it to the equilibrium surface mass balance (SMB) distribution obtained from a glacier flow model simulation (Zwinger et al., 2014) of the SBB valley, using the identical surface DEM. Assuming a non-varying glacier surface, s, the equilibrium SMB is given by where (u ice , v ice , w ice ) T are the Cartesian components of the surface ice velocity vector. The SMB distribution obtained by Zwinger et al. (2014) is depicted in Fig. 7. The sensitivity of the wind-impact distribution with respect to the inflow profile was assessed by rerunning with the slightly altered inflow conditions, as presented in Eqs. (6) and (7). We then assess the sensitivity of the solution to bottom surface topography by using the smoothed DEM. Finally, we assess the effects of the LGM topography on the wind fields.

Variations of inflow profiles
The first simulation is conducted using the nose-shaped profile given by the cubic polynomial (Eq. 5). We run the simulation for slightly longer than 20 min of physical time. Figure 8 shows the velocity magnitude at 5 m above ground as well as on a vertical cross section cutting through the area above the northern inner BIA from east to west, displaying snapshots starting at 500 s (when the first gusts reach the valley bottom) in steps of 250 to 1250 s (the end of the simulation). An animated GIF (showing steps of 50 s) of this as well as of all runs conducted can be found in the Supplement of this article. From these snapshots it is clear that fast flow is focused on the inner BIAs. The wind speeds are consistently high above the ablation area (which largely coincides with the inner BIA) found from the glacier flow simulation, indicated by the yellow polygon in Fig. 7. In particular the close correspondence in the position of the high-velocity streak with the northern boundary of the larger inner BIA at 1250 s into the simulation is striking. Figure 9 shows the average as well as the standard deviation of the velocity for the time window from 950 to 1050 s of the simulation, with peak values for both properties concentrated over the inner northern BIA. The latter can be interpreted as a measure of velocity fluctuations on the scale of the finest mesh size (vertically in the range of metres), which in LES represent the eddies containing the highest turbulent kinetic energy (Axelsen and van Dop, 2009). Figure 10 shows results from similar simulations with altered inflow profiles as in Eqs. (6) and (7) at 1000 s into the simulation. Although the patterns obtained with these pro-    Fig. 6 with the results in Figs. 8 and 10 we can clearly see that the closer the maximum flow speed of the inflow profile is to the surface, the greater the focusing of the fast flow towards the BIA seems to be. From this we conclude that the high wind speed close to the terrain creates the effect of redirection.

Variations in the smoothness of the terrain
We are interested in how small-scale topographic features may affect the flow details and hence the impact of the DEM accuracy. To that end we smoothed the original DEM using the previously described algorithm and applied the noseshape inflow defined by Eq. (5). The left panel in Fig. 11 (and animated GIF in the Supplement) shows the velocity magnitude 5 m above ground as well as on a vertical cross section from east to west at 1000 s into the simulation. One striking difference from all runs conducted on the original terrain is that at 1000 s into the simulation the katabatic front has advanced about 1 km further into the valley. The general im-pression -even more supported by the series in the animated GIF than by the single snapshot -is that the areas of highimpact wind velocities are less focused on the inner BIAs than the reference run (Fig. 8).

LGM simulation
The surface elevation inside the valley of SBB may have been considerably higher (Näslund et al., 2000;Hättestrand and Johansen, 2005) at the LGM. We increased the ice thickness by a constant 200 m inside and 100 m outside the valley and smoothed the surface to resemble glaciation of present-day nunataks (i.e. ice-free areas on the bounding mountains of SBB valley). The right panel in Fig. 11 shows the distribution of the velocity magnitude 5 m above ground and on a vertical cross section from east to west. In contrast with all other simulations conducted on the present-day DEM -either in its original or smoothed form -the high-velocity winds in the LGM simulation are distributed more randomly over the valley and not focused over present-day BIAs. In fact, the highest velocities at 1000 s into the simulation are located at the north-western end of the valley, close to the position of the outer BIA.

Conclusions
By conducting local wind simulations above the Scharffenbergbotnen valley we show that katabatic surface jets coming from the east produce high-impact wind speeds at exactly the location of the inner blue ice area. These types of wind were probably directly experienced in the field during Jan-uary 2007 when a highly localized event caused a clearing of surface firn, likely to be several years old, from underlying bare ice. During this event, simultaneously taken direct and AWS observations support the dominant role of local katabatic events rather than synoptic storm events in formation of the Scharffenbergbotnen BIA.
Sensitivity analysis to changes in vertical wind-profile shapes and DEMs revealed that the topography is mainly responsible for the distribution of the high-wind-impact features. For all simulation runs with the present-day DEM, the derived locations of both fast wind flow and recirculation areas match well with the location of the inner BIA in the Scharffenbergbotnen valley. In contrast, a simulation conducted with increased ice thickness imposed on the presentday DEM (resembling conditions at Late Glacial Maximum) revealed that snow and firn removal focused at the position of the inner BIA most likely occurred only after the lowering of the ice in the SBB valley. Only then did the present-day shape of the surrounding rocks surface from beneath the retreating and lowering ice sheet. In other words, the inner BIA at Scharffenbergbotnen is younger than the LGM. This is in agreement with the findings of Zwinger et al. (2014), who came to similar conclusion by looking at age/depth horizons of a full-stress ice flow model.
We suggest that other BIAs in similar mountain regimes can be effectively generated by katabatic jets and that these jets are produced independently of synoptic storms. It is conceivable that other BIAs may be controlled by extremely powerful synoptic storms that occur from time to time and scour snow and firn from the surface of the BIAs.
The Supplement related to this article is available online at doi:10.5194/tc-9-1415-2015-supplement.