Introduction

Owing to their low albedo (0.5–0.6) compared with that of dry snow (0.8–0.9) (Reference Bintanja and van den BroekeBintanja and Van den Broeke, 1995; Reference RasmusRasmus, 2006), blue-ice areas (i.e. areas of negative mass balance) absorb approximately twice as much solar energy as is absorbed by snow (Reference BintanjaBintanja, 1999, Reference Bintanja2000). At elevations above 1000 m a.s.l, with lower air temperature, blue-ice areas are wind-induced and rarely contain internal water. However, in low-elevation blue-ice areas (LEBIAs) extensive subsurface ponds can form (Reference Winther, Elvehøy, Bøggild, Sand and ListonWinther and others, 1996). Because of the low heat conductivity of ice (2.2 W m⁻¹ K⁻¹ for pure ice at 0°C), which is still large compared with that of snow, the absorbed radiative energy escapes very slowly, forming a subsurface temperature maximum. This is called the solid-state greenhouse effect and is described by Reference ColbeckColbeck (1989) and Reference Brandt and WarrenBrandt and Warren (1993). If the temperature reaches 0°C, the ice can begin to melt. The subsurface meltwater ponds formed subsequently contain a large amount of heat that could cause melting in the surrounding snow cover if released into it. Reference Liston and WintherListon and Winther (2005) calculated the meltwater flux from blue-ice areas and found it to be 57.5 km³ a⁻¹. They concluded that the meltwater does not influence the mass balance of the ice sheet because it probably refreezes locally. Observations of meltwater flow in blue-ice areas have been made by Reference Winther, Elvehøy, Bøggild, Sand and ListonWinther and others (1996) and Reference PhillipsPhillips (1998) indicating that the water does move around. Reference Winther, Jespersen and ListonWinther and others (2001), using satellite mosaic imagery of Antarctica, have found the present spatial extent of blue-ice areas to be ∼60000 km² each for LEBIAs and wind-induced blue-ice areas.

Variations in the spatial extent of blue-ice areas have been used previously as indicators of climate change (e.g. Reference Orheim and LucchittaOrheim and Lucchitta, 1990), but the occurrence of subsurface ponds has not. Their occurrence, spatial extent and depth could be used as indicators of warming (Reference Rasmus and BeckmannRasmus and Beckmann, 2007), but the processes involved have not been adequately observed or modelled. The colour of the blue-ice area changes when liquid water is present, because of the decreased scattering from water compared with ice. This makes the concept of remotely sensing the internal melting feasible to some extent.

Previous models of subsurface melting in blue-ice areas (Reference Liston, Winther, Bruland, Elvehøy and SandListon and others, 1999a,Reference Liston, Bruland, Winther, Elvehøy and Sandb) have focused on the thermodynamic aspects of the phenomenon in one vertical dimension, which implies the necessary assumption of a semi-infinite homogeneous blue-ice area. This assumption was relaxed in the study on the hydrodynamics and evolution of subsurface ponds by Reference Rasmus and BeckmannRasmus and Beckmann (2007), which used a model with a two-dimensional (2-D) set-up. We refer to this 2-D model as the RB model. Rasmus and Beckmann found that the movement of water transports heat downward in the pond. Vertical convection and a weak overturning circulation are generated, which acts to homogenize the fluid and transport warmer water downward, thereby causing additional melting at the base of the pond.

In the Reference Rasmus and BeckmannRasmus and Beckmann (2007) study, the ice began to disintegrate internally during a 50 year integration with an air-temperature trend of 3°C (100 a)⁻¹. The RB model had a domain that was only 4.8 m deep and it used a no-flux lower boundary condition. This was also the case in the model of Reference Liston, Winther, Bruland, Elvehøy and SandListon and others (1999a,Reference Liston, Bruland, Winther, Elvehøy and Sandb). The use of a no-flux lower boundary condition assumes the domain is deep enough to damp out temperature variations. If this is not the case, then heat will start to build up because it is not able to exit the model domain. The disintegration of the ice found by Reference Rasmus and BeckmannRasmus and Beckmann (2007) could therefore have been due to heat not being able to escape downwards, so that a lower boundary condition that would let the heat through could lead to a different result.

The RB model also used a spectrally constant bulk attenuation coefficient for ice and water, together with a spectrally constant albedo at the ice surface. The constant attenuation coefficients could also have caused an overestimation of the radiative heating rates as indicated by Reference Brandt and WarrenBrandt and Warren (1993) and by Reference Warren, Brandt, Grenfell and McKayWarren and others (2002) in their study of sea-ice formation in snowball Earth conditions. Finally, because the RB model is 2-D, it is possible to study the effects of spatial variations in surface properties on this phenomenon. Blue-ice areas can have large spatial variations of surface properties (Fig. 1). The effect of these on the absorption of radiation and the flow of heat is shown schematically in Figure 2.

Fig. 1. LEBIA below Basen nunatak on which the Finnish Antarctic research station, Aboa, is located (73°03′ S, 13°25′ W). Plogen nunatak is in the background. Notice the spatial variations in surface reflectance visible in the image. Drinking water is being pumped from the blue-ice area. Photograph taken by K. Rasmus, January 2000, during the Finnish Antarctic Research Programme (FINNARP) 99 expedition to Antarctica.

Fig. 2. Schematic representation of ice partially covered by a thin layer of snow. More radiation comes to the middle part of the ice and it heats up more than the sides, which are shaded by the high-albedo surface. This leads to a net flow of heat from the middle to the sides.

Our purpose here is to study the effects of the different boundary conditions and the attenuation coefficient used in a thermo-hydrodynamic model on the subsurface melting in an idealized LEBIA. The intention is not to simulate any existing blue-ice area.

Our objectives are:

1. 1. To test the changes made to the RB model. These include a radiating lower boundary condition in the model together with a spectral albedo and spectral attenuation coefficients.

2. 2. To quantify the amount of subsurface meltwater available in the ice during the austral summer after 50 years with and without a specified trend in the air temperature.

3. 3. To study the effects of spatial variations in surface albedo on the subsurface melt-pool formation.

This study is a further development of the RB model and it starts by describing the changes made to the model. Results from the different model calculations with the model and the integrations with and without the trend in air temperature are then presented. The paper then concludes with a discussion on the results and their significance for the Antarctic ice sheet.

The Model

The model has been introduced in detail by Reference Rasmus and BeckmannRasmus and Beckmann (2007). It consists of three parts: (1) an optical part that describes the penetration of shortwave radiation into the ice/water mixture; (2) a thermodynamic part that describes the conductive heat transfer in ice and water; and (3) a hydrodynamic part that describes the water–ice matrix circulation within the subsurface pond. The changes made to the model and the model configuration are described in more detail here, together with a more detailed description of the optical part.

In this study the no-flux lower boundary condition used in the RB model is changed to a radiating boundary condition of the form:

(1)

where the phase velocity, c, is taken from the previous gridpoints in the z direction by assuming that the flux was the same at the bottom boundary as at the previous gridcell. Equation (1) states that no heat was allowed back into the domain from below, which assumes a semi-infinite ice sheet. No other changes were made to the thermodynamic part.

The optical part of the model was changed to incorporate multiband optical properties. An eight-band optical model was used with the following bands: 200–400, 400–600, 600–700, 700–800, 800–1000, 1000–1200, 1200–1600 and 1600–3000 nm. Here, we use the word ‘spectral’ as a synonym of ‘multiband’.

The optical simulation started by calculating the global irradiance, Q ₀, at each time-step. This depends mainly on the cloud cover, which was set to a constant 0.5, and the solar elevation. Downwelling spectral irradiance, E, is then calculated by multiplying Q ₀ with the normalized spectral distribution, Q _s(λ). To calculate the amount of radiation penetrating into the ice, the downwelling irradiance is multiplied by (1−α), where α is the albedo. A spectral albedo for ice was used to calculate the penetrating spectral irradiance; otherwise a bulk value was used for total irradiance. Figure 3 shows Q _s(λ) and the spectral albedo for blue ice. The spectral albedo for blue ice is from Reference RasmusRasmus (2006). A bulk albedo value of 0.85 is used for snow.

Fig. 3. Spectral albedo for blue ice and the normalized spectral distribution, Q _s, of incoming solar radiation.

Starting from the definition of the total albedo in which the spectral albedo is integrated over the whole solar spectrum,

(2)

and with the knowledge that the downwelling irradiance is E _d(λ) = Q ₀ Q _s(λ), the total albedo is seen to be

(3)

When the spectral albedo in Figure 3 is integrated in this way, the bulk total albedo is found to be ∼0.50.

The irradiance is then transferred through the ice using the exponential decay law for irradiance, which in a general form is

(4)

where K(λ) is the spectral attenuation coefficient and λ is the wavelength. For water the irradiance is taken from values for the clearest oceanic waters presented by Reference Rasmus, Granéli and WängbergRasmus and others (2004). The spectral attenuation coefficient for blue ice is a combination of blue-ice measurements presented by Reference Rasmus, Granberg, Kanto, Kärkäs, Lavoie and LeppärantaRasmus and others (2003) and model calculations of Reference Liston, Bruland, Winther, Elvehøy and SandListon and others (1999b). The spectral attenuation coefficients for blue ice and water are shown in Figure 4. The Mie model was used by Reference Liston, Winther, Bruland, Elvehøy and SandListon and others (1999a) to calculate the optical properties dynamically, but this approach is not used in this study because the Mie model requires grain size as an input variable and the grain size of clear ice can be ambiguous. A discussion of this issue can be found in Reference Warren, Brandt, Grenfell and McKayWarren and others (2002). For the purposes of this idealized study, one spectrum of attenuation coefficient for ice is sufficient.

Fig. 4. Spectral attenuation coefficients for blue ice and water.

Solar radiation in the medium was split into a part through ice and a part through water within each gridcell (Fig. 5). In both cases, exponential decay with depth was assumed. The water and ice were assumed to be in parallel vertical columns within each gridcell.

Fig. 5. Representation of the gridcell used in the model. The cell is divided vertically into a water part and an ice part. The radiation decays in the two parts separately.

The radiative contribution to the heat flux was then computed using the exponential decay presented in Equation (4) separately for ice and water through each gridcell. The attenuation coefficients were constant with depth, which simplified the equation. The results were then summed at the bottom of each gridcell to achieve the total radiation at the bottom:

(5)

where A is the ice concentration, that is A = 1 means 100% ice and A = 0 means 100% liquid water. K _i and K _w are the attenuation coefficients for ice and water, respectively.

No modifications were made to the hydrodynamic part of the model except for the relaxation of the hydrodynamically active domain boundary, which is the ice concentration at which the hydrodynamics of the water become active in the model. The ice concentration is the ice to water volume ratio. The domain boundary was changed from an ice concentration of 0.0 to 0.2, because it was assumed that at these ice concentrations the ice is present in such small quantities that it does not affect the internal friction (viscosity) of the water. The boundary is checked at each time-step and the hydrodynamical equations are evaluated when necessary.

The surface heat budget is computed taking into account the shortwave incoming radiation (including diurnal and seasonal variations), albedo, longwave incoming and out-going radiation as well as sensible and latent heat fluxes, following the parameterization proposed by Reference ZillmanZillmann (1972). The forcings are the same as in the original RB model. Because this is an idealized process study, the model uses highly schematic surface forcings for temperature, wind, humidity and cloud cover to capture the typical seasonal cycle in Antarctica at 70° S. These are described in more detail by Reference Rasmus and BeckmannRasmus and Beckmann (2007). The forcings used are simplifications that do not take into account the stratification of the atmospheric boundary layer, an issue that needs to be addressed in future versions of the model.

Air temperature (T _a) and dewpoint temperature (T _d) have a functional dependence:

(6)

(7)

where t is the time in days from the beginning of the first year. In Equations (6) and (7), the second term on the right is the yearly temperature cycle and the third term is the temperature trend of 3°C (100 a)⁻¹ which was applied to some of the 50 year integrations.

The model configuration features a 2-D (x-z) domain with a 9.6 m lateral extent with 100 gridpoints, and a 9.6 m vertical extent with 300 gridpoints. This was twice the extent of the Reference Rasmus and BeckmannRasmus and Beckmann (2007) model domain. The spatial resolution was 3.2 cm in the vertical and 6.4 cm in the lateral direction. The domain is considered periodic in the lateral x direction and a radiating boundary condition was used at the bottom for heat. The size of the domain was dictated by total calculation time after the choice of gridcell size, which came from the demand for numerical stability with a sensible time-step, was made. The time-step used was 20 s.

Model Calculations

The albedo was set to the blue-ice albedo in the middle part of the domain. Initially, the width of the area was 8.9 m. The albedo was set to the snow albedo (0.85) at the edges, mimicking a thin layer of snow on the surface of the ice (Fig. 2). The initial conditions corresponded to a −10°C isothermal glacier temperature and no liquid water.

The model years consisted of 12 months with 30 days in each. The integrations started at the beginning of the austral summer in October of the first year. The different model calculation set-ups are summarized in Table 1.

Table 1. Parameters for each model run. The hydrodynamically active domain was at an ice concentration of 0.2 for all

To isolate the effects of the bottom boundary condition, the experiment of Reference Rasmus and BeckmannRasmus and Beckmann (2007) was repeated with the same bulk albedos and attenuation coefficients. First, the model was run without a trend in air temperature and then with a trend of 3°C (100 a)⁻¹. These are runs 1a and 1b in Table 1. Second, the bulk attenuation coefficient was changed to a spectral attenuation coefficient and three runs were made, first with two bulk albedos of 0.5 and 0.65, and then with the spectral albedo shown in Figure 3. These are runs 2a, 2b and 2c in Table 1.

The model was then run with a spatially varying albedo to study the effect of spatial surface variations on the subsurface melting phenomenon. The centre had a spectral blue-ice albedo and the edges had a bulk snow albedo. When a part of the domain is in shade, the amount of absorbed radiation within the ice becomes smaller. This leads to a heat flow from the parts that are not in shade (Fig. 2). The width of the area with the blue-ice spectral albedo was set to values between 0.0 and 3.8 m and the albedo outside this area was set to the snow albedo (0.85). This is run 3 in Table 1.

For the final experiment the model was run with a spatially and spectrally constant bulk surface albedo that varied from 0.2 to 1.0, and a spectral attenuation coefficient. This is run 4 in Table 1.

Results

The midday bottom temperatures for the first model integrations (runs 1a and 1b), testing the lower boundary conditions, are shown in Figure 6 with both temperature scenarios, that is constant mean air temperature and with the trend. The results show that, without the trend in the air temperature, it takes 7 years for the model to stabilize at the bottom depth. Maximum temperatures are obtained at the end of August, showing that the heatwave from the surface propagates in 5 years to the 9.6 m depth. This depth is still not enough to damp out the heatwave completely, but because the temperature variations remain relatively stable for the no-trend air temperature case, it would seem that the bottom boundary condition is working. There is no significant leakage or build-up of heat at the bottom depth.

Fig. 6. Midday bottom (9.6 m) temperatures for both air-temperature scenarios for the 50 year integration. See runs 1a and 1b in Table 1.

The melting pond maximum-extent results in runs 1a and 1b for the middle of the model domain (Fig. 7) show that the melting is in the top 100 cm of the ice and that no disintegration of the ice pond occurs within the 50 years that the model was run even with the trend in air temperature. The effect of the trend is seen as a slight deepening of the pond with time, and a slight shallowing of the solid ice layer at the surface. At the surface, the model stabilizes after 1 year.

Fig. 7. Time evolution of the top and bottom limits of the melt pool at its maximum extent for both scenarios, with and without a trend in air temperature, for the 50 year integration. See runs 1a and 1b in Table 1. A 1 year evolution figure of ice concentration can be seen in Reference Rasmus and BeckmannRasmus and Beckmann (2007).

Changing the albedo from bulk values to spectral values, but still using the spectral attenuation coefficient, has quite a large effect on the subsurface melting (runs 2a, 2b and 2c). The situation in the middle of the domain at the end of March of the second year (Fig. 8) shows that the minimum ice concentrations for the cases with a spectral albedo and with the bulk albedo value of 0.65 are very similar, being just over 0.2. The overall depth of the pond is larger (nearly 3 m) for the bulk albedo case. The depth of the pond in the spectral albedo (0.65) case is 0.5 m less. However, when the corresponding bulk albedo value of 0.50 is used, the pond is deepened considerably relative to the spectral albedo, and the minimum ice concentration value goes to 0.0. Because the hydrodynamically active domain is set to an ice concentration of 0.2, the motion of the water has affected the deepening in this case. The surface layer containing solid ice is the same thickness in all cases. It has been kept frozen by radiative cooling, sublimation and heat loss due to the sensible heat flux.

Fig. 8. Ice-concentration profiles from the middle of the domain at the time of minimum ice concentration with two bulk albedos (0.50 and 0.65) and a spectral albedo. See runs 2a, 2b and 2c in Table 1. Notice the dramatic decrease in subsurface melting when the bulk albedo of 0.5 is changed to a spectral albedo.

The values of ice concentration from 0.5 m depth for bulk albedos varying between 0.20 and 1.00 are shown in Figure 9 (run 3). The albedo dictates the amount of radiation available for penetration into the ice. Because the penetration of radiation is described in a spectral way, the net result of this albedo increase is to decrease the amount of shortwave radiation penetrating into the ice per waveband. The results show that the threshold for subsurface melting is 0.70–0.80. It is not possible for snow to have extensive internal melting because it has an albedo above 0.80. When the surface albedo is below 0.35, the subsurface pool persists through the austral winter. At 0.5 m depth an ice concentration below 0.1 is only possible with an albedo <0.50. As the albedo decreases from 1.00, the timing of the minimum ice concentration moves forward. The range in time of total freeze-up is not very large within an albedo value of 0.40–0.70. These calculations used a spectral attenuation coefficient but they used bulk values for the albedo.

Fig. 9. Time series of ice concentration from 0.5 m depth in the model domain. See run 3 in Table 1. The first 2 years of the model run are shown.

The effect of the width of the blue-ice area on the ice concentration (run 4) at 0.5 m depth is only seen up to a width of 0.64 m (10 gridpoints) after which the ice concentration reaches a minimum value of 0.23 (Fig. 10). The edges of the area with the blue-ice albedo cause a slight increase in ice concentration up to a distance of ∼0.5 m from the edge. This is because heat can escape into the ice shaded by the thin layer of snow (mimicked by the surface having a higher albedo).

Fig. 10. Ice concentration values from 0.5 m depth fordifferent widths of the area with the blue-ice albedo. See run 4 in Table 1. The widths in the legend are in metres.

Discussion

In their study of subsurface melting in LEBIAs, Reference Rasmus and BeckmannRasmus and Beckmann (2007) found that even though the surface remains well below freezing and the top 20 cm remains as solid ice, the heat that has reached >2 m depth cannot escape during the winter and slowly leads to a total subsurface disintegration of the ice with the prescribed air-temperature increase. In their case, the subsurface pond acted as a thermal barrier for the heat below because the water was not able to transport heat up to the surface. In this study with the radiating lower boundary condition, the pond thickness increased in the global change scenario but at a much lower rate (Fig. 7) compared with the original study with the RB model. The majority of the increase was due to deepening. Even after 50 years with an air temperature increase of 1.5°C, the liquid water did not reach the surface, and the disintegration of the LEBIA did not occur. This could be due to the heat being able to pass through the bottom boundary or in part due to the almost twice as large domain compared with the original RB model. However, the domain size is probably not a major factor because the radiation penetration and absorption is mainly a one-dimensional (1-D) process and, therefore, the same in all gridpoints.

The results show that the use of a radiating boundary can decrease the build-up of heat and excessive melting. This radiation of heat through the bottom boundary is illustrated by the variation in bottom temperature that follows the variation in annual air temperature (Fig. 6). If the LEBIA is taken to be semi-infinite in the vertical direction, the conclusion of this study would be that the LEBIA will not disintegrate after 50 years even with the specified increase in air temperature. The heat would just continue down into the ice infinitely. Because the ice sheet is not semi-infinite the fate of the heat transported through the bottom boundary should be considered. Is it transported laterally, in which case it could re-enter the LEBIA at a later stage? Is it trapped further down in the ice, in which case the ice could still disintegrate after some time or at a different location? It would then be important to add external lateral transport to the model. Even though the timescale of the Reference Rasmus and BeckmannRasmus and Beckmann (2007) break-up may not be accurate, the break-up of the blue-ice area should still occur because the capacity of the glacier to absorb heat is finite and the heat that flows through the bottom boundary must end up somewhere.

Reference Winther, Jespersen and ListonWinther and others (2001) found the total spatial extent of LEBIAs to be approximately 60 000 km². The results in this study show that after 50 years the depth of the water layer within the ice was 1 m during the summer for the constant climate scenario, increasing to 1.4 m for the global change scenario (Fig. 7). In the figures, the limits for the water layer are set at an ice concentration of 0.5. These results mean that during the summer in Antarctica, LEBIAs have the potential to form a maximum of 60 km³ of meltwater which would increase to 84 km³ if the air temperature were to rise by the specified amount of 3°C (100 a)⁻¹. Because the melt pools do not extend to the edges of LEBIAs, the real volumes of water are less than these. The formation of subsurface water in a LEBIA and its discharge into the surrounding snow will have an impact on the ice sheet because of the large latent heat content of the water. Even though blue-ice areas constitute only 1% of the area of Antarctica and only half of that is melt-induced, this phenomenon can still be important because it is a mechanism by which mass is moved around.

A bulk albedo cannot be substituted for the spectral albedo without causing a difference in the subsurface melting. A spectral albedo corresponding to a bulk albedo of 0.5 causes almost the same melting as a bulk albedo of 0.65. If the bulk value is used for each spectral band of incoming solar radiation, the amount of penetrating shortwave radiation in the visible region is much higher than in the spectral case. A surface with an albedo of 0.5 absorbs 1.4 times as much radiation as a surface with an albedo of 0.65. This means that the spectral distribution of radiation together with the spectral attenuation coefficient plays an important role in this process.

When comparing runs with bulk optical coefficients (runs 1a and 1b; Fig. 7) with runs with spectral values (runs 2b and 2c; Fig. 8), it is seen that the solid surface extends deeper in run 2b, but also the vertical thickness of the area of meltwater is deeper in the spectral case. This is because the shortwave radiation has been split up into wavebands and each waveband has been treated with a different attenuation coefficient. There is less infrared radiation in the spectrum, and this is the part with the largest attenuation coefficient. The visible radiation, with a smaller attenuation coefficient, is able to penetrate deeper into the ice. This explains the greater depth of the upper surface of the water layer when using a spectral attenuation coefficient.

The variation of albedo in blue ice is 0.5–0.7 (Reference RasmusRasmus, 2006). The results show that at these albedo values there is always the potential for subsurface melting in a LEBIA. However, even a thin snow cover can increase the albedo to a value above 0.7 making it impossible for substantial subsurface melting to occur. On the other hand, when the snow begins to thin above a blue-ice area, the subsurface melting can start even before the snow cover is totally gone. Within the range of naturally occurring albedos, the subsurface meltwater pool cannot persist over the austral summer.

Snow and its spatial distribution have a large effect on subsurface melting. First, by having a much higher albedo it is able to diminish the melting directly below it and, if it is spatially very heterogeneous, it can shade the blue ice next to it and diminish the melting further. The spatial distribution of the snow on top of blue-ice areas needs to be taken into account for more accurate meltwater flux calculations. In addition to snow, the surface of the ice can be composed of grains with distinct grain boundaries and it can contain air bubbles. Both of these increase the albedo, but they also increase the attenuation coefficient (e.g. Reference Warren, Brandt, Grenfell and McKayWarren and others, 2002). This has not been taken into account in this study, so with a snow surface albedo the radiation is able to penetrate deeper than in a real snow-covered case.

The results shown here for the snow shading calculations (Fig. 10) are valid for a directionally uniform diffuse incident flux (Fig. 2). In a real situation where the incident flux consists of both a direct (arriving at an oblique angle) and a diffuse component, the shading effect of the snow cover on the blue-ice area can be even greater than presented here, especially at low solar elevations. The resulting length scale needs to be taken into account (e.g. when calculating meltwater fluxes from blue-ice areas obtained from satellite images) because it is much smaller than the resolution of the images currently available with remote sensing. The layer of water within the ice changes the spectrum of radiation scattered to the sensor. This can be used to track melting with a satellite remote-sensing instrument.

The ice concentrations obtained in this study in all the modelled cases are qualitatively similar to those obtained by previous modelling (Reference Liston, Winther, Bruland, Elvehøy and SandListon and others, 1999a,Reference Liston, Bruland, Winther, Elvehøy and Sandb; Reference Rasmus and BeckmannRasmus and Beckmann, 2007). The previous models without a hydrodynamical component, however, were only able to produce a minimum ice concentration of 80%. As expected, the initial melting and liquid water formation is mainly a 1-D (vertical) process. This means that horizontal transport of water and impurities (e.g. dust, sand or soot) in the subsurface water layer must be due to other processes (e.g. larger-scale horizontal advection). These processes are beyond the scope of this study but are necessary components of future studies.

Many other scenarios could be investigated: changing temperatures might be accompanied by changes in surface wind, cloudiness and humidity. The local surface topography (flat vs sloping), the presence or absence of cracks and crevasses, as well as the thickness of snow will determine the specific evolution of a given LEBIA. Also, a more detailed surface heat budget based on measured values is necessary to calibrate and validate the model results to some specific LEBIAs.

A proper validation with real data of the model is important before it is made more complicated. Low-elevation blue-ice data in Antarctica are still scarce, however, but automatic weather stations from the vicinity of blue-ice areas exist and they are able to produce real and usable forcings for the model.

The model, and especially the optical part of it, still requires further development. Implementation of a two-stream multiband approximation for the penetration of the solar radiation can increase the accuracy of the solar radiation profile, especially close to the top surface. The addition of an external lateral flow would enable the model to calculate the spreading of substances. For climate-change studies, the air-temperature increase could be changed to an increase in longwave radiation from the atmosphere.

Conclusions

We have made sensitivity studies for a model capable of calculating subsurface melting in a LEBIA. Thermodynamical and hydrodynamical aspects in a model study of subsurface ponds in LEBIAs have been taken into account before, but this study shows what the effects are of varying the boundary conditions and optical properties. The conclusions of this study are:

1. 1. The lower boundary has a large effect on the subsurface melting results. A radiating boundary condition diminishes the melting. It was found that a bulk albedo cannot be substituted for the spectral albedo without causing a difference in the subsurface melting. A spectral albedo corresponding to a bulk albedo of 0.50 causes almost the same melting as a bulk albedo of 0.65. This causes more radiation to be available for penetrating into the ice.

2. 2. When spectral radiation was passed through the medium with a spectral attenuation coefficient, the pool of water was thicker and the top surface was deeper than when using bulk optical properties.

3. 3. The meltwater pool was 1 m thick after a 50 year integration. This increased to 1.5 m when the trend in air temperature was applied. The limits of the pond were set at an ice concentration of 0.5. For a typical atmospheric warming scenario of 1.5°C in 50 years, subsurface changes were found. The surface remains frozen and largely unchanged during the 50 year integration. The ice did not disintegrate even after a 50 year integration.

4. 4. The bulk albedo threshold for subsurface melting to occur is 0.7–0.8. The subsurface ponds were deep enough to survive the summer when the bulk albedo was below 0.4. In this set-up, the hydrodynamics play a minor role at large albedo values but have an important role when the albedo reaches the threshold value. The corresponding integrated spectral albedos will be lower than these.

5. 5. Shading by snow has a direct effect on subsurface melting when the snow areas are <0.6 m apart. This effect is increased by negligible melting below the area covered by snow. These points have remote-sensing implications. Because this study uses a downwelling flux arriving from the zenith, the length scale is probably less in nature.

More work is still required and a real blue-ice case would be interesting to study for validation purposes.