The Cryosphere Interaction between ice sheet dynamics and subglacial lake circulation : a coupled modelling approach

Subglacial lakes in Antarctica influence to a large extent the flow of the ice sheet. In this study we use an idealised lake geometry to study this impact. We employ a) an improved three-dimensional full-Stokes ice flow model with a nonlinear rheology, b) a three-dimensional fluid dynamics model with eddy diffusion to simulate the basal mass balance at the lake-ice interface, and c) a newly developed coupler to exchange boundary conditions between the two individual models. Different boundary conditions are applied over grounded ice and floating ice. This results in significantly increased temperatures within the ice on top of the lake, compared to ice at the same depth outside the lake area. Basal melting of the ice sheet increases this lateral temperature gradient. Upstream the ice flow converges towards the lake and accelerates by about 10% whenever basal melting at the ice-lake boundary is present. Above and downstream of the lake, where the ice flow diverges, a velocity decrease of about 10% is simulated.


Introduction
During the last decades our knowledge on subglacial lake systems has greatly increased.Since the discovery of the largest subglacial lake, Lake Vostok (Oswald and Robin, 1973;Robin et al., 1977), more than 270 other lakes have been identified so far (Siegert et al., 2005;Carter et al., 2007;Bell, 2008;Smith et al., 2009).Plenty of efforts have been undertaken to reveal partial secrets subglacial lakes may hold, mostly referring to Lake Vostok.For instance the speculation that extremophiles in subglacial lakes may be encountered (e.g., Duxbury et al., 2001;Siegert et al., 2003) has Correspondence to: M. Thoma (malte.thoma@awi.de)been nurtured by microorganisms discovered in ice core samples (Karl et al., 1999;Lavire et al., 2006).These samples originate from the at least 200 m thick accreted ice, drilled at the Russian research station Vostok (Jouzel et al., 1999).Discussions about the origin and history of the lake (Duxbury et al., 2001;Siegert, 2004;Pattyn, 2004;Siegert, 2005) are still ongoing.It is still unknown whether Subglacial Lake Vostok existed prior to the Antarctic glaciation and if it could have survived glaciation, or whether it formed subglacially after the onset of glaciation.The impact of subglacial lakes on the flow of the overlying ice sheet has been analysed (Kwok et al., 2000;Tikku et al., 2004, e.g.,) and modelled (Mayer and Siegert, 2000;Pattyn, 2003;Pattyn et al., 2004;Pattyn, 2008), but self-consistent numerical models, which include the modelling of the basal mass balance, are lacking at present.Several numerical estimates and models (West and Carmack, 2000;Williams, 2001;Walsh, 2002;Mayer et al., 2003;Thoma et al., 2007Thoma et al., , 2008b;;Filina et al., 2008) as well as laboratory analogues (Wells and Wettlaufer, 2008) of water flow within the lake have been carried out, permitting insights into the water circulation, energy budget and basal processes at the lake-ice interface.Finally, observations (Bell et al., 2002;Tikku et al., 2004) and numerical modelling (Thoma et al., 2008a) of accreted ice at the lakeice interface gave insights about its thickness and distribution over subglacial lakes.
Early research assumed that subglacial lakes were isolated systems, while only recently Gray et al. (2005), Wingham et al. (2006) andFricker et al. (2007) found evidence that Antarctic subglacial lakes are connected, hence forming an extensive subglacial hydrological network.Several plans to unlock subglacial lakes exist (Siegert et al., 2004;Inman, 2005;Siegert et al., 2007;Schiermeier, 2008;Woodward et al., 2009) and valuable knowledge will be available as soon as direct samples of subglacial water and sediments are taken.
Published by Copernicus Publications on behalf of the European Geosciences Union.
However, current knowledge about the interaction between subglacial lakes and the overlying ice sheet is lacking.The most important parameter exchanged between ice and water is heat.The exchange of latent heat associated with melting and freezing dominates heat conduction, but the latter process has also be accounted for in subglacial lake modelling (Thoma et al., 2008b).Melting and freezing is closely related to the ice draft which varies spatially over subglacial lakes (Siegert et al., 2000;Studinger et al., 2004;Tikku et al., 2005;Thoma et al., 2007Thoma et al., , 2008aThoma et al., , 2009a)).The ice draft, and hence the water circulation within the lake, is maintained by the ice flow across lakes.Without this flow, lake surfaces would even out (Lewis and Perkin, 1986).On the other hand, a spatially varying melting/freezing pattern will have an impact on the overlying ice sheet as well.In order to get an insight in the complex interaction processes between the Antarctic ice sheet and subglacial lakes, this study for the first time couples a numerical ice-flow and a lake-flow model with a simple asynchronous time-stepping scheme to overcome the problem of different adjustment time scales.
We describe the applied ice-flow model RIMBAY and lakeflow model ROMBAX in the first two Sections, respectively.Each section starts with a general description of the particular model, describes the applied boundary conditions, and the results of the uncoupled model runs.The newly developed RIMBAY-ROMBAX-coupler RIROCO is introduced in Sect.4, where we also discuss the impact of a coupled icelake system on the ice flow and lake geometry.

General description
The ice sheet model RIMBAY (Revised Ice sheet Model Based on frAnk pattYn) is based on the work of Pattyn (2003), Pattyn et al. (2004) and Pattyn (2008).Within this thermomechanically coupled, three-dimensional, full-Stokes ice model, a subglacial lake is represented numerically by a vanishing bottom-friction coefficient of β 2 = 0, while high friction in grounded areas is represented by a large coefficient; in this study we use β 2 = 10 6 .The choice of the latter parameter has no influence on the model results as long as it is much larger than zero.
The constitutive equation, governing the creep of polycrystalline ice and relating the deviatoric stresses τ to the strain rates ε, is given by Glen's flow law: ε = Aτ n (e.g., Pattyn, 2003), with a temperature dependent rate factor A = A(T ).Here we apply the so-called Hooke's rate factor (Hooke, 1981).Experimental values of the exponent n in Glen's flow law vary from 1.5 to 4.2 with a mean of about 3 (Weertman, 1973;Paterson, 1994); ice models traditionally assume n = 3.However, a simplified viscous rheology with n = 1 stabilises and accelerates the convergence behaviour of the implemented numerical solvers.Therefore, previous subglacial lake simulations within ice models were limited to viscous rheologies (Pattyn, 2008).

Model setup and boundary conditions
The model domain used in this study is a slightly enlarged version of the one presented in Pattyn (2003): It consists of a rectangular 168 100 km 2 large domain with a model resolution of 5 km (resolutions of 2.5 km and 10 km are also used for comparison) and 41 terrain-following vertical layers.The surface of the initially 4000 m thick ice sheet has an initial slope of 2% (similar to Lake Vostok, Tikku et al., 2004) from left (upstream) to right (downstream).An idealized circular lake with a radius of about 48 km and an area of about 7200 km 2 is located in the center of the domain, where a 1000 m deep cavity modulates the otherwise smoothly sloped bedrock (similar to Pattyn, 2008).The lake's maximum water depth of about 600 m, resulting in a volume of about 1840 km 3 (the lake's geometry is also indicated in Figs. 3, 4).
We apply a constant surface temperature of −50 • C at the ice's surface, a typical value for central Antarctica (Comiso, 2000).The bottom layer boundary temperature depends on the basal condition: above bedrock a Neumann boundary condition is applied, based on a geothermal heat flux of 54 mW/m 2 (a value suggested for the Lake Vostok region by Maule et al., 2005).Assuming isostasy above the subglacial lake, the bottom layer temperature is at the pressuredepended freezing point.This temperature, prescribing a Dirichlet boundary condition, depends on the local ice sheet thickness (T b = −H ×8.7×10 −4 • C/m ,e.g., Paterson, 1994).Accumulation and basal melting/freezing are ignored during the initial experiments in this Sect., where no coupling is applied.In the subsequent coupling-experiments, basal melting and freezing at the lake's interface, as modelled by the lake flow model, is accounted for (Sect.4).
The lateral boundary conditions are periodic, hence values of ice thicknesses, velocities, and stresses are copied from the upstream (left) side to the downstream (right) side of the model domain and vice versa.The same applies to the lateral borders along the flow.The initial integration starts from an isothermal state.The basal melt rate is set to zero over bedrock as well as over the lake.

Model improvements
Compared to Pattyn (2008) we improved the model RIMBAY in two significant ways to allow a numerically stable and fast representation of the more realistic non-linear flow law with an exponent of n = 3: First, a gradual increase of the friction coefficient β 2 at the lake's boundaries is considered (Fig. 1a).Physically, this smoothing can be interpreted as a lubrication of the ice sheet base on the grounded side and as stiffening due to debris on the lake side of the lake's edge.Our experiments have shown, that a slight prescribed β 2 -smoothing The Cryosphere, 4, 1-12, 2010 www.the-cryosphere.net/4/1/2010/coefficient of (1/0) is suited to decrease the integration time significantly and stabilises the numerical results.In this notation, (1) indicates one lubricated node on the bedrock side of the boundary and (0) indicates no stiffed (debris) node on the lake side.If the number of lubricated nodes is reduced to zero, the model's integration time increases without a significant impact on the velocity field.If the number of stiffed nodes over the lake is increased, lower velocities over the lake are achieved.However, for higher model resolutions than used in this study, higher β-values should be considered to make the transition more realistic.
Second, we implement a three-dimensional Gaussian-type filter to smooth the viscosity as well as the vertical resistive stress with a variable filter width.Without this filter, the numerics becomes unstable.Figure 1b shows a one-dimensional equivalent to the implemented three dimensional filter.Figure 2 compares three different smoothing parameter results on the friction coefficient β 2 as well as the Gaussian-type filter impact on the viscosity for three different filter widths.Sensitivity experiments have shown, that with a slight transient β 2 smoothing at the lake edges combined with a gentle (quadratic) Gaussian filter, as shown in Fig. 2b, numerical stability is achieved.

Results
The standard experiment is a thermomechanically coupled full-Stokes (FS) model with a horizontal resolution of 5 km.
A quasi-steady state is reached after 300 000 years. Figure 3a shows the lake's position and depth in the center of the model domain.The frictionless boundary condition over the lake results in an increased velocity, not only above the lake, but also in it's vicinity (Fig. 3b).From mass conservation it follows, that the (vertically averaged) horizontal velocity converges towards the lake from upstream and diverges downstream.The flattened ice sheet surface over the lake is a consequence of the isostatic adjustment.This effect is also visible in the geometry of the profiles in Fig. 4 and in particular in Fig. 5.The inclined lake-ice interface is maintained by the constant ice flow over the lake and hence, independent of a possible basal mass exchange.
Figure 3b shows the vertically averaged horizontal velocity and Fig. 4a a vertical cross section of the velocity along the flow at the lake's center at y = 200 km.Over the lake, the ice sheet behaves like an ice shelf, featuring a vertically constant velocity.Towards the lake, the surface velocity increases by more than 70% from about 0.7 m/a to about 1.2 m/a.The largest velocity gradients occur in the vicinity of the grounding line (Figs. 3b,4a,5).Convergence of ice towards the lake results in an accelerated ice flow which undulates at the lake's edges because of significant changes in the ice thickness and ice sheet surface gradients (yellow and red line in Fig. 5, respectively).The surface velocity (green line in Fig. 5) accelerates towards the lake until the ice sheet surface flattens.The velocity increase in deeper layers close to the bedrock (blue line in Fig. 5) is suspended just before the frictionless lake interface is reached, because of a strong increase of the ice thickness filling the trough.The maximum velocity is reached across the lake (Figs. 3b,4a,5), where the basal friction is zero.The vertical temperature profile (Fig. 4b) is nearly linear, as accumulation and basal melting are neglected.Geothermal heat flux is not sufficient to melt the bottom of the grounded ice, but over the lake the freezing point is maintained by the boundary condition.This results in submerging isotherms.At the downstream grounding line a slight overshoot of the upwelling isotherm is observed.

Robustness of the results
To investigate the impact of the horizontal resolution on the model results, simulations with a coarser 10 km as well as a finer 2.5 km grid resolution have been performed.Within the coarser resolution most features are reproduced well, but there is a significant impact on ice velocities, in particular along the grounding line where the differences reaches about 10%.In general, the model with the higher spatial resolution has increased velocities along the flowlines across the lake and decreased velocities outside.In addition, the local The Cryosphere, 4, 1-12, 2010 www.the-cryosphere.net/4/1/2010/velocity maximum at the grounding line (Fig. 5) cannot be resolved with the 10 km resolution.The finer 2.5 km model resolution needs a much longer integration time, without producing significant differences in the results compared to the intermediate 5 km grid resolution.
We also investigated, whether a higher order model, neglecting resistive stress and vertical derivatives of the vertical velocity (see Pattyn 2003, Saito et al. 2003, Marshall 2005  3 Lake flow model "ROMBAX"

General description
To simulate the water flow in the prescribed subglacial lake we apply "ROMBAX", a terrain-following, primitive equations, three-dimensional, fluid dynamics model (e.g., Griffies, 2004)."ROMBAX" simulates the interaction between ice and subjacent water in terms of melting and freezing, according to heat and salinity conservation and the pressure dependent freezing point at the interface (Holland and Jenkins, 1999).The model uses spherical coordinates and has been applied successfully to ice-shelf cavities (e.g., Grosfeld et al., 1997;Williams et al., 2001;Thoma et al., 2006) as well as to subglacial lakes (Williams, 2001;Thoma et al., 2007Thoma et al., , 2008a 6a).Across the following gradient from abou side of ice flow to about 24 m results from the modelled tem (Fig. 4b) and is used as input

Results
The to reproduce the results obdel.All other aspects of the etrization are kept identical.erate impact on the ice flow: subglacial lake are about 5% kes model.In the vicinity of -Stokes terms decreased, but ines.Although the FS-model in the numercial code it conder model, hence all further ull-Stokes model.

X"
n the prescribed subglacial further details).The horizontal resolution (0.025 • ×0.0125 • , about 0.7 × 1.4 km), the number of vertical layers ( 16), as well as the horizontal and vertical eddy diffusivities (5 m 2 /s and 0.025 cm 2 /s, respectively) are adopted from a model of subglacial Lake Concordia (Thoma et al., 2009a).In a model domain of about 170 × 88 × 16 grid cells the circulation within the lake as well as the melting and freezing rates at the lake-ice interface are calculated.At the bottom of the lake a geothermal heat flux of 54 mW/m 2 , consistent with the ice-flow model's boundary condition, is applied.Previous subglacial lake simulations of Lake Vostok (Thoma et al., 2007(Thoma et al., , 2008a;;Filina et al., 2008), Lake Concordia (Thoma et al., 2009a), or Lake Ellsworth (Woodward et al., 2009) used a prescribed average heat conduction into the ice (Q Ice = dT /dz × 2.1 W/(K m)), based on borehole temperature measurements and thickness temperature-gradient estimates.The availability of the results of the thermomechan- geometry as well as the parametrization are kept identical.Our experiments indicate a moderate impact on the ice flow: calculated velocities across the subglacial lake are about 5% lower compared to the full-Stokes model.In the vicinity of the lake, the impact of the full-Stokes terms decreased, but is still enhanced along the flowlines.Although the FS-model implements more mathematics in the numercial code it converges faster than the higher order model, hence all further studies are performed with the full-Stokes model.

General description
To simulate the water flow in the prescribed subglacial lake we apply ROMBAX, a terrain-following, primitive equations, three-dimensional, fluid dynamics model (e.g., Griffies, 2004).ROMBAX simulates the interaction between ice and subjacent water in terms of melting and freezing, according to heat and salinity conservation and the pressure dependent freezing point at the interface (Holland and Jenkins, 1999).The model uses spherical coordinates and has been applied successfully to ice-shelf cavities (e.g., Grosfeld et al., 1997;Williams et al., 2001;Thoma et al., 2006) as well as to subglacial lakes (Williams, 2001;Thoma et al., 2007Thoma et al., , 2008a,b;,b;Filina et al., 2008;Thoma et al., 2009a,b,c;Woodward et al., 2009)  only parameter exchange beween both models was the unilateral initilisation of the lake flow model ROMBAX with the 370 modelled geometry and temperature gradient of the ice flow model RIMBAY (indicated by the gray lines in Figure 7).The real coupling procedure starts, when the results of ROMBAX are reinserted into RIMBAY (indicated by the black lines in Figure 7).

Results
From ROMBAX the basal mass balance at the ice sheet-lake interface is considered for subsequent initialisations of RIMBAY.Melting dominates freezing (which is negligible) and hence the ice sheet is loosing mass.Note that the molten ice does 380 not affect the lake's volume, neither in the ice-sheet model nor in the lake-flow model, as this is constant per definition.This mass imbalance can be interpreted as a virtual constant water flow out of the lake without any feedback to the ice sheet.The coupler RIROCO applies restarts from 385 previous model runs.Consequently the models reach their new quasi-steady state after a significant shorter integration time.Ice sheet volume in the model domain is decreasing by about 600 km 3 per 100 000 years, equivalent to about 3.5 m thickness.Ice draft reduction in the lake-flow model within 390 100 000 years is shown in Figure 8a.Most ice is lost in the center of the lake (up to 8 m), and the area of maximum mass loss is slightly shifted to the upstream side of the lake.Consequently, the water column thickness increases where most ice is molten and decreases where less ice is molten (Fig- 395 ure 8b).However, the ice-thickness reduction and slope adjustment is too small to change the surface pressure on the lake water, and hence the water flow, significantly.Just a few iteration cycles are needed to bring this coupled lake-water system into a quasi-steady state.

400
Several impacts of bottom melting on ice on top of the lake are observed: First, the temperature gradient at the ice sheet's bottom is increased.This results in an increase of heat conduction into the ice by about 22% (Figures 8d), compared to a model run without basal melting (Figures 6a).Consequently, 405 more heat is extracted from the lake and the modelled average melting decreases by about 7% (Figures 8c).Second, ice sheet thickness above the lake is reduced.This increases the surface gradient towards the upstream part and decreases the surface gradient towards the downstream part.Hence, the ice 410 flow upstream accelerates and decelerates downstream (Figure 9a).The magnitude of the ice velocity change is about 10%.Third, bottom melting removes mass and a vertical downward velocity follows from mass conservation.This advects colder ice from the surface towards the bottom, re-415 sulting in a relative cooling of up to -2.1 • C (Figure 9b) above the lake.This negative temperature (compared to the model

Model setup and boundary conditions
The bedrock topography and the ice draft, needed for the lake-flow model ROMBAX, is obtained from the output (Sect.2.4) of the ice-flow model RIMBAY (see Sect. 4 for further details).The horizontal resolution (0.025 • ×0.0125 • , about 0.7 × 1.4 km), the number of vertical layers (16), as well as the horizontal and vertical eddy diffusivities (5 m 2 /s and 0.025 cm 2 /s, respectively) are adopted from a model of subglacial Lake Concordia (Thoma et al., 2009a).In a model domain of about 170 × 88 × 16 grid cells the circulation within the lake as well as the melting and freezing rates at the lake-ice interface are calculated.At the bottom of the lake a geothermal heat flux of 54 mW/m 2 , consistent with the ice-flow model's boundary condition, is applied.Previous subglacial lake simulations of Lake Vostok (Thoma et al., 2007(Thoma et al., , 2008a;;Filina et al., 2008), Lake Concordia (Thoma et al., 2009a), or Lake Ellsworth (Woodward et al., 2009) used a prescribed average heat conduction into the ice (Q Ice = dT /dz × 2.1 W/(K m)), based on borehole temperature measurements and thickness temperature-gradient estimates.The availability of the results of the thermomechanical ice-sheet model RIMBAY permits a spatially varying Q Ice (Fig. 6a).Across the lake's center, a general draft-following gradient from about 27 mW/m 2 on the upstream side of ice flow to about 24 mW/m 2 on the downstream side results from the modelled temperature distribution in the ice (Fig. 4b) and is used as input for the lake-flow model.

Results
The initial model run starts with a lake at rest.After about 200 years a quasi-steady state is reached.The circulation within the lake is shown in Fig. 6b-c.The vertically integrated mass transport stream function (with a strength of about 1.3 mSv, 1 mSv=1000 m 3 /s) as well as the zonal overturning (about 0.3 mSv) show a two-gyre structure, while the meridional overturning (about 1.3 mSv) indicates just one anticyclonic gyre.The strength of the mass transport is between those modelled for Lake Vostok and those for Lake Concordia (Thoma et al., 2009a), and hence reasonable for subglacial lakes.There is only a slight ice draft slope from about 4006 m to 3927 m across the 90 km of the lake (Figs. 5,4).This results in a decrease of melting from about 12 mm/a in the West to a negligible freezing along the eastern shoreline (Fig. 6d).A significant amount of freezing would only be modelled if the ice draft would have a steeper slope.

General description
The bedrock topography and the ice draft, necessary to set up the geometry for the lake-flow model ROMBAX (Sect.3.2), was obtained from the modelled output geometry of the iceflow model RIMBAY (Sect.2.4).In addition, the thermodynamic boundary condition of the heat conduction into the ice was calculated from the temperature gradient at the ice sheet's bottom.The left part of Fig. 7 (gray lines) shows schematically the performed operations.A coordinate conversion is necessary as RIMBAY uses Cartesian coordinates while ROMBAX is based on spherical coordinates.
The modelled melting and freezing rates (Fig. 6d) are considered to replace the previously zero-constrained lower boundary condition in the ice-flow model RIMBAY (indicated by the central triangle within the yellow area in Fig. 7).Again a coordinate transformation is necessary.To speed up the integration time, ice geometry, ice flow, and ice temperature from the initial model run are reused (indicated by the black-lined triangle within the blue area within Fig. 7).An additional process has to be considered since the lake's area may change during an ice-model run.As the basal mass balance (melting/freezing) depends on the former lake-model results and cannot be calculated during a specific model run, extrapolation of neighbouring values may be necessary (indicated by the embedded yellow oval in the blue area of Fig. 7).
In each coupling cycle, successive initialisations of the lake-flow model ROMBAX are performed with the slightly changed ice draft and water column thickness (indicated by the right triangle in the yellow area of Fig. 7), as well as the temperature field of the previous model-run (indicated by the central triangle in the orange area of Fig. 7).Because the lake flow model does not permit dynamically changing geometry, temperatures of emerging nodes have to be extrapolated from neighbouring nodes (indicated by the embedded yellow oval in the orange area of Fig. 7).It is not necessary to reuse (and possibly extrapolate) the water circulation, as this value is based on the tracer (density) distribution and converges quickly.
The coupling mechanism, including initial starts of the individual models, parameter exchanges, coordinate transformations, and restart-initialisations are embedded into and controlled by the RIMBAY-ROMBAX-Coupler RIROCO.
In Sect.2.4 and Sect.3.3 the results of the initial runs of the individual models have been described.So far the only parameter exchange beween both models was the unilateral initilisation of the lake flow model ROMBAX with the modelled geometry and temperature gradient of the ice flow model RIMBAY (indicated by the gray lines in Fig. 7).The real coupling procedure starts, when the results of ROMBAX are reinserted into RIMBAY (indicated by the black lines in Fig. 7).

Results
From ROMBAX the basal mass balance at the ice sheet-lake interface is considered for subsequent initialisations of RIM-BAY.Melting dominates freezing (which is negligible) and hence the ice sheet is loosing mass.Note that the molten ice does not affect the lake's volume, neither in the ice-sheet model nor in the lake-flow model, as this is constant per definition.This mass imbalance can be interpreted as a virtual constant water flow out of the lake without any feedback to the ice sheet.The coupler RIROCO applies restarts from previous model runs.Consequently the models reach their new quasi-steady state after a significant shorter integration time.Ice sheet volume in the model domain is decreasing by about 600 km 3 per 100 000 years, equivalent to about 3.5 m thickness.Ice draft reduction in the lake-flow model within 100 000 years is shown in Fig. 8a.Most ice is lost in the center of the lake (up to 8 m), and the area of maximum mass loss is slightly shifted to the upstream side of the lake.Consequently, the water column thickness increases where most ice is molten and decreases where less ice is molten (Fig. 8b).However, the ice-thickness reduction and slope adjustment is too small to change the surface pressure on the lake water, and hence the water flow, significantly.Just a few iteration cycles are needed to bring this coupled lake-water system into a quasi-steady state.
Several impacts of bottom melting on ice on top of the lake are observed: First, the temperature gradient at the ice sheet's bottom is increased.This results in an increase of heat conduction into the ice by about 22% (Fig. 8d), compared to a model run without basal melting (Fig. 6a).Consequently, more heat is extracted from the lake and the modelled average melting decreases by about 7% (Fig. 8c).Second, ice sheet thickness above the lake is reduced.This increases the surface gradient towards the upstream part and decreases the surface gradient towards the downstream part.Hence, the ice flow upstream accelerates and decelerates downstream (Fig. 9a).The magnitude of the ice velocity change is about 10%.Third, bottom melting removes mass and a vertical downward velocity follows from mass conservation.This advects colder ice from the surface towards the bottom, resulting in a relative cooling of up to −2.1 • C (Fig. 9b) above the lake.This negative temperature (compared to the model run without melting, Fig. 4b) is advected downstream by the horizontal velocity.The lake-ice interface itself is cooled by less than 0.007 • C, as it is maintained at the pressure-dependent freezing point.The cooling, visible at the lateral boundaries in Fig. 9b, result from the applied periodic boundary condition.It is an artifact of the numerical representation of the setup and not discussed here.

Summary
Observations from space show that large subglacial lakes have a significant impact on the shape and dynamics of the Antarctic Ice Sheet as they flatten the ice sheets surface (Siegert et al., 2000;Kwok et al., 2000Kwok et al., , 2004;;Leonard et al., 2004;Tikku et al., 2004;Siegert, 2005).These regions depict a change in surface slope due to the isostatic adjustment of the ice sheet.This, combined with the lacking bottom friction, results in an observable redirection of ice flow.In this study, we apply a newly coupled ice sheet-lake flow model on an idealized ice sheet-lake configuration to investigate the feedbacks between the individual systems.The full-Stokes ice sheet model (Pattyn, 2008) is improved to handle a more realistic non-linear rheology.The lake-flow model is based on a three dimensional fluid dynamics model with eddy diffusivity and simulates the lake flow as well as the mass balance at the lake-ice interface.This model, previously only applied to lake-exclusive studies with a prescribed ice thickness and bedrock, now receives its geometry directly from the ice flow model.In addition, heat flux from the lake into the ice sheet is an exchange parameter.Besides other external forcing fields, the ice-flow model receives the basal mass balance at the interface from the dynamic lake-flow model.
In order to stabilise the ice-flow model numerically with a nonlinear rheology, a slight Gaussian filtering of the ice The ice flow converges where an ice-shelf type flow servation requires a downs of the deceleration when th Melting at the lake-ice inte eration upstream and decre as well as downstream, bec This impacts on the veloci and downstream of the la corresponds to about 20 tim temperature at the ice sheet pressure melting point, and based ice in the vicinity.A ing at the lake-ice interface the lake.Advection transpo stream, and hence the tem have impacts on the ice flo Our idealised configurati an important impact of the on the ice sheet's dynamic coupled model will focus o Lake Vostok and its glacia vestigate the flow dynamic with observations.tion between ice sheet and subglacial lakes 9 The ice flow converges and accelerates towards the lake, where an ice-shelf type flow structure establishes.Mass conservation requires a downstream divergence of flow, because of the deceleration when the flow reaches the bedrock again.Melting at the lake-ice interface increases the ice flow acceleration upstream and decreases the ice flow on top of the lake as well as downstream, because of the reduced surface slope.This impacts on the velocity can be traced about 75 km upand downstream of the lake in this specific setting, which corresponds to about 20 times the ice sheet thickness.The temperature at the ice sheet bottom on top of the lake is at the pressure melting point, and hence warmer than the bedrockbased ice in the vicinity.According to our simulation, melting at the lake-ice interface reduces the temperature on top of the lake.Advection transports the relatively colder ice downstream, and hence the temperature dependent rheology will

Summary
Observations from space show that large subglacial lakes have a significant impact on the shape and dynamics of the Antarctic Ice Sheet as they flatten the ice sheets surface (Siegert et al., 2000;Kwok et al., 2000;Leonard et al., 2004;Tikku et al., 2004;Siegert, 2005).These regions depict a change in surface slope due to the isostatic adjustment of the ice sheet.This, combined with the lacking bottom friction, results in an observable redirection of ice flow.In this study, we apply a newly coupled ice sheet-lake flow model on an idealized ice sheet-lake configuration to investigate the feedbacks between the individual systems.The full-Stokes ice sheet model (Pattyn, 2008) is improved to handle a more realistic non-linear rheology.The lake-flow model is based on a three dimensional fluid dynamics model with eddy diffusivity and simulates the lake flow as well as the mass balance at the lake-ice interface.This model, previously only applied to lake-exclusive studies with a prescribed ice thickness and bedrock, now receives its geometry directly from the ice flow model.In addition, heat flux from the lake into the ice sheet is an exchange parameter.Besides other external forcing fields, the ice-flow model receives the basal mass balance at the interface from the dynamic lake-flow model.
In order to stabilise the ice-flow model numerically with a nonlinear rheology, a slight Gaussian filtering of the ice viscosity is necessary.
The ice flow converges and accelerates towards the lake, where an ice-shelf type flow structure establishes.Mass conservation requires a downstream divergence of flow, because www.the-cryosphere.net/4/1/2010/The Cryosphere, 4, 1-12, 2010 of the deceleration when the flow reaches the bedrock again.
Melting at the lake-ice interface increases the ice flow acceleration upstream and decreases the ice flow on top of the lake as well as downstream, because of the reduced surface slope.This impacts on the velocity can be traced about 75 km upand downstream of the lake in this specific setting, which corresponds to about 20 times the ice sheet thickness.The temperature at the ice sheet bottom on top of the lake is at the pressure melting point, and hence warmer than the bedrockbased ice in the vicinity.According to our simulation, melting at the lake-ice interface reduces the temperature on top of the lake.Advection transports the relatively colder ice downstream, and hence the temperature dependent rheology will have impacts on the ice flow beyond the lake.
Our idealised configuration of an ice-lake system indicates an important impact of the interaction between both systems on the ice sheet's dynamics.Next applications of this new coupled model will focus on realistic configurations, such as Lake Vostok and its glacial drainage system.In order to investigate the flow dynamical impact which can be compared with observations.

Fig. 2 .Fig. 2 .
Fig. 1. a) Friction coefficient β 2 (scaled to the maximum of 10 6 ) along the central x-axis for an unsmoothed and two test-smoothed cases.The smoothing coefficient (x/y) represents the number of nodes smoothed outside and inside of the lake's border, respectively.b) Weight factors for a Gaussian-type filter with different width from one (red) to five (cyan) depending on the distance to the central node located at zero.

Fig. 1 .Fig. 1 .Fig. 2 .
Fig. 1.(a) Friction coefficient β 2 (scaled to the maximum of 10 6 ) along the central x-axis for an unsmoothed and two test-smoothed cases.The smoothing coefficient (x/y) represents the number of nodes smoothed outside and inside of the lake's border, respectively.(b) Weight factors for a Gaussian-type filter with different width from one (red) to five (cyan) depending on the distance to the central node located at zero.

Fig. 2 .
Fig. 2. Logarithm of viscosity η of the surface layer and bedrock-friction parameter β 2 for an idealized model.(a) No β 2 -smoothing or Gaussian filtering.(b) Slight β 2 -smoothing (1/0) and Gaussian filtering with a filter width of two.(c) Strong β 2 -smoothing (3/2) and Gaussian filtering with a filter width of five.Note that the Figs.represent the initial conditions for a FS-model with the applied filter and β 2 -smoothing.

Fig. 3 .
Fig. 3. Results for a full-Stokes ice dynamic model with a horizontal resolution of 5 km.(a) Lake depth (color), ice sheet surface elevation (dashed contours), and ice-flow velocity (black arrows).(b) Vertically averaged horizontal velocity of the standard full-Stokes experiment.
,b;Filina et al., 2008; Thoma et al., 2009a,b,c;Woodward et al., 2009).3.2Model setup and boundary conditionsfurther details).The horizonta about 0.7 × 1.4 km), the num well as the horizontal and ver and 0.025 cm 2 /s, respectively of subglacial Lake Concordia model domain of about 170 × tion within the lake as well as at the lake-ice interface are ca lake a geothermal heat flux o the ice-flow model's boundar vious subglacial lake simulat etal., 2007, 2008a; Filina e  (Thoma et al., 2009a), or Lak 2009) used a prescribed avera (Q Ice = dT /dz × 2.1 W/(K m) ture measurements and thickn mates.The availability of the ical ice-sheet model "RIMBA Q Ice (Figure initial model run starts w 200 years a quasi-steady stat within the lake is shown in tegrated mass transport stream about 1.3 mSv, 1 mSv=1000 m turning (about 0.3 mSv) show M. Thoma et al.: Modelled interaction between ice sheet and subglacial lakes Ice ).(b) Vertically integrated mass transport stream function, positive values indicate a clockwise lockwise circulation.(c) Zonal (from west to east) and meridional (from south to north) overturning.(d)

Fig. 6 .
Fig. 6.(a) Heat conduction into the ice (Q Ice ).(b) Vertically integrated mass transport stream function, positive values indicate a clockwise circulation, negative values an anticlockwise circulation.(c) Zonal (from west to east) and meridional (from south to north) overturning.(d) Basal mass balance.

Fig. 7 .
Fig. 7. Couple scheme schematics.The upper part represents the ice-flow model RIMBAY, the lower part the lake-flow model ROMBAX.In the middle the parameters exchanged by the coupler RIROCO are shown.The gray lines in the left part of the figure represent the initial start-up sequence (the two initial model runs), dashed loops indicate cycles (successive model runs) that may repeat an apriori unspecified number of times.Two additional yellow ovals indicate the individual model needs regarding successive restarts/coupling: The ice flow calculated by RIMBAY may result in new lake nodes.For these nodes the basal mass balance is calculated by averaging adjacent nodes.For ROMBAX a modified geometry desires the extrapolation of temperatures (and other tracers) from a previous model run to skip the time-consuming spin-up process. 375

Fig. 7 .
Fig. 7. Couple scheme schematics.The upper part represents the ice-flow model RIMBAY, the lower part the lake-flow model ROMBAX.In the middle the parameters exchanged by the coupler RIROCO are shown.The gray lines in the left part of the Fig. represent the initial start-up sequence (the two initial model runs), dashed loops indicate cycles (successive model runs) that may repeat an apriori unspecified number of times.Two additional yellow ovals indicate the individual model needs regarding successive restarts/coupling: The ice flow calculated by RIMBAY may result in new lake nodes.For these nodes the basal mass balance is calculated by averaging adjacent nodes.For ROMBAX a modified geometry desires the extrapolation of temperatures (and other tracers) from a previous model run to skip the time-consuming spin-up process.

Fig. 8 .
Fig. 8. Impact of melting after 100 000 years on (a) the ice draft, (b) the water column thickness, and the geometry difference between the lake-flow model restart and the initial geometry.(d) Indicates the (a) the ice draft, (b) the water column thickness, and (c) the basal mass balance.Shown is ake-flow model restart and the initial geometry.(d) Indicates the heat conduction into the ice (Q Ice ).rical representation of the at large subglacial lakes hape and dynamics of the en the ice sheets surface 000, 2004;Leonard et al.,  005).These regions depict he isostatic adjustment of h the lacking bottom fricection of ice flow.In this ice sheet-lake flow model guration to investigate the

Fig. 8 .
Fig. 8. Impact of melting after 100 000 years on (a) the ice draft, (b) the water column thickness, and (c) the basal mass balance.Shown is the geometry difference between the lake-flow model restart and the initial geometry.(d) Indicates the heat conduction into the ice (Q Ice ).