Thermo-hydro-mechanical processes in fractured rock formations during a glacial advance

The paper examines the coupled thermo-hydromechanical (THM) processes that develop in a fractured rock region within a fluid-saturated rock mass due to loads imposed by an advancing glacier. This scenario needs to be examined in order to assess the suitability of potential sites for the location of deep geologic repositories for the storage of high-level nuclear waste. The THM processes are examined using a computational multiphysics approach that takes into account thermo-poroelasticity of the intact geological formation and the presence of a system of sessile but hydraulically interacting fractures (fracture zones). The modelling considers coupled thermo-hydro-mechanical effects in both the intact rock and the fracture zones due to contact normal stresses and fluid pressure at the base of the advancing glacier. Computational modelling provides an assessment of the role of fractures in modifying the pore pressure generation within the entire rock mass.


Introduction
The longevity constraints with regard to the safety of deep geologic sequestration of high level radioactive waste suggest that the conventional scientific approaches to the investigations that involve laboratory and field studies must be complemented by approaches that will allow for predictions on timescales that are beyond conventional geological engineering activities involving underground facilities (Laughton et al., 1986;Chapman and McKinley, 1987;Selvadurai and Nguyen, 1997;Rutqvist et al., 2005;Alonso et al., 2005). Current concepts for deep geologic storage require an assessment of the geological setting containing a repository to account for geomorphological processes that can occur over timescales of several thousands of years. The major geomorphological process that is identified over this timescale is glaciation. Our attention is restricted to a geologic setting that has been investigated in connection with the DECOVALEX project, and the domain of interest incorporates a system of fractures that are hydraulically connected but mechanically sessile under the influence of the glaciation loading. The domain of interest, which contains the set of fractures, is derived from the international DECOVALEX III project (Chan and Stanchell, 2008). The highly idealized geosphere model is based on data from the Whiteshell Research Area (WRA) in Manitoba, Canada. The studies by Stanchell (2005, 2008) provide the database of material properties in the simulations conducted on stationary aspects of glacial loads with specific emphasis on the hydro-mechanical modelling that uses the continental-scale model of the Laurentide Ice Sheet developed by Boulton et al. (2004). The dominant fracture network closely resembles the simplified version of an actual fracture network found in a small (approx. 10 km × 10 km) subregion of the Canadian Shield at WRA.
The purpose of this paper is to examine coupled thermohydro-mechanical (THM) problems for a fluid-saturated rock mass that contains fracture zones and is subjected to advancing glacial loads. The mechanical response of a fractured rock mass to a glacial cycle was studied in the papers by Selvadurai and Nguyen (1995, 1999 and Steffen et al. (2014); it was shown that fracture zones can lose stability or slip during or after the deglaciation period of the cycle. The glacial cycle scenario was also included in the DE-COVALEX studies

of a fractured rock
Published by Copernicus Publications on behalf of the European Geosciences Union. 2168 A. P. S. Selvadurai et al.: Thermo-hydro-mechanical processes in fractured rock formations mass and the current study extends this work to include THM effects and a more accurate treatment of the glacial loading. A secondary objective of the study is to assess the capabilities of a computational multiphysics finite element code COMSOL Multiphysics® in examining a model domain that consists of an intact rock mass and fracture zones, both of which can exhibit THM processes.
The formulation of a coupled hydro-mechanical (HM) problem for fully saturated geological materials was presented by Biot (1941) and reformulated by several investigators including Rice and Cleary (1976) and Hamiel et al. (2005). The HM coupling gives rise to the Mandel-Cryer effect (Selvadurai and Yue, 1994;Selvadurai and Shirazi, 2004), which can explain the momentary increase in the pressure head in an aquifer when the ground water is drained upwards along a fault (Stanislavsky and Garven, 2003) and can predict changes in the fluid pressure induced by the slip of geological faults (Hamiel et al., 2005). Thermohydro-mechanical effects have also been of interest to many geoenvironmental processes including the geologic disposal of heat-emitting nuclear waste, frictional heating of faults (Lachenburch, 1980;Rice, 2006), geothermal energy extraction (Dickson and Fanelli, 1995), and ground freezing resulting from buried chilled-gas pipelines (Selvadurai et al., 1999a, b). Solutions to pure fluid flow in heterogeneous formations and coupled thermo-poroelastic problems for fluidsaturated geological materials have been obtained by several researchers (Booker and Savvidou, 1985;McTigue, 1986;Nguyen, 1995, 1997;Nguyen and Selvadurai, 1995;Selvadurai, 1996;Khalili and Selvadurai, 2003;Selvadurai and Selvadurai, 2010;Suvorov, 2012, 2014). Recent mathematical and computational studies of THM behaviour of fluid-saturated media with both elastic and elasto-plastic skeletal behaviour are given by Suvorov (2012, 2014) and experimental manifestations of the thermo-poroelastic Mandel-Cryer effect are also given by Najari and .
The behaviour of a rock mass can be influenced by the presence of fractures and faults. Coupled THM behaviour of fractured rocks has been studied by several investigators and extensive references to this work can be found in areas related to geosciences and geomechanics (Noorishad et al., 1984;Selvadurai and Nguyen, 1995;Nguyen et al., 2005;Rutqvist et al., 2002;Guvanasen and Chan, 2000;Chan and Stanchell, 2008;Tsang et al., 2009). Thermohydro-mechanical processes in fractured rock formations can be analysed using two approaches: in the first approach, by modelling discrete fractures and specifying their locations within a host rock, which is modelled as a fracturefree medium; and in the second approach by introducing the influence of fractures implicitly through the derived overall constitutive equations for the fractured medium.
Within a discrete fracture modelling approach, three distinct finite element formulations can be identified: special interface or joint elements Nguyen, 1995, 1999;Ng and Small, 1997;Nguyen and Selvadurai, 1998;Guvanasen and Chan, 2000;Steffen et al., 2014), the embedded manifold approach (Guvanasen and Chan, 2000;Juanes et al., 2002;Graf and Therrien, 2008;Erhel et al., 2009), and the conventional or direct approach, in which the fractures are modelled with the finite elements of the same spatial order; e.g. 3-D fractures are modelled with 3-D finite elements (Stanislavsky and Garven, 2003;Sykes et al., 2011). The multiphysics finite element code COMSOL Multiphysics® used in this study has an efficient 3-D mesh generator, well suited for discretizing complex geometries containing distinct narrow regions, such as fractures. Thus, in this paper, the THM problems investigated in connection with glacial loading will be examined using the conventional or direct approach.
In this study the rock mass is assumed to be an isotropic thermo-poroelastic domain with a network of dominant fractures that is integral with the surrounding intact rock. The intact rock is assumed to contain small-scale joints, minor fractures, pores and voids, and thus they are not explicitly included into the model. The dominant fracture network represents large-scale faults and closely resembles the simplified version of the real fracture network found in a small (approx. 10 km × 10 km) subregion of the Canadian Shield described in the NWMO (Nuclear Waste Management Organization) report by Chan and Stanchell (2008). The fracture network is placed at the interior of the rock mass domain and remote from its boundaries. Consequently, two distinct regions can be identified in the given rock mass: an interior region with fractures and an exterior region void of fractures. This allows for the observation of the THM response of the fractured region, particularly as the glacier advances, without any dominant influences of the boundaries.
In this work attention is restricted to a glacier of very large length which is initially at rest and then starts to move along the bedrock surface. The glacier motion or advance is assumed to be very small compared to the glacier length. Since the glacier loading is highly non-uniform near the glacier front and very flat near its centre, the major change in the response of the underlying bedrock caused by a small glacier advance is expected to occur near the initial location of the glacier front. Therefore, we will place the rock mass model domain containing fractures near the glacier front initial position and observe the response of the rock and fractures there. The size of the model domain is also chosen to be very small compared to the glacier length although larger sizes would naturally provide for better accuracy of the solution. Since the nonlinear effects are not included in the present study, the initial fields existing in the rock before motion of the glacier has started can conveniently be set to zero, and then only the change in the response of rock caused by the glacier advance can be observed. The tectonic stress effects of glacial isostatic adjustment (GIA) (see e.g. Peltier, 2004) are not accounted for in the current modelling, although it is an important consideration in a more complete investigation of the glacial loading problem. The crustal movements can continue for several millennia due to creep and viscoelastic effects (Walcott, 1970(Walcott, , 1976McNutt and Menard, 1978;Selvadurai, 1979) and a glaciation episode needs to be considered in the context of GIA that needs to be accurately estimated to provide a reference state.
The paper is organized as follows. Mathematical description of the thermo-poroelasticity problem is given in Sect. 2 whereas the finite element model is described in detail in Sect. 3. A few numerical tests that validate the model are presented in Sect. 4 and Sect. 5 contains the main computational results that illustrate the influence of glacial loading and temperature change on the development of fluid pressure, velocity, temperature, displacement and mean effective stress in the entire rock mass and within the individual fractures.

Constitutive models
For a linear elastic isotropic fluid-saturated rock, the total stress tensor σ ij is related to the infinitesimal strain tensor ε ij , the fluid pressure p, and the temperature T via the constitutive equation: In Eq. (1), ε V is the volumetric strain, K D and G D denote respectively the bulk and shear modulus of the drained rock, α s is the thermal expansion coefficient of the solid phase of the rock, and α is the Biot coefficient defined as α = 1−K D /K s , where K s is the bulk modulus of the solid phase. The effective stress tensor σ ij is defined by and it represents the stress supported by the solid skeleton.
For an isotropic fluid-saturated porous medium, Darcy's law can be written: where v i is the spatially averaged fluid velocity, k is the permeability of the rock, and µ is the dynamic viscosity. Also, ρ f denotes the fluid density and g is the gravitational acceleration; here it is assumed that the gravitational acceleration vector is pointing in the negative direction of the vertical z axis. The rocks encountered in deep geological settings are relatively impervious in their intact state and the permeability parameter is an important input to the modelling exercise. The estimation of permeability in such low permeability rocks is, however, a non-routine exercise. Transient tests (Selvadurai and Carnaffan, 1997;Selvadurai and Jenner, 2012;Selvadurai and Najari, 2013;Najari and Selvadu-rai, 2014;Selvadurai et al., , 2011) are used to estimate the fluid transport properties of low permeability rocks and even in these tests, factors such as presence of air in the pressurized fluid cavity and in the accessible pore space of the rock, the degree of saturation, residual pore fluid pressures and the stress state -induced micro-cracks and damage -can influence the estimated permeability (Selvadurai, 2004(Selvadurai, , 2009aSelvadurai and Głowacki, 2008;Selvadurai and Ichikawa, 2013;Selvadurai, 2012, 2014;Najari, 2013, 2015). We assume that the heat transfer in the system is through conduction only and the fluid flow velocity both in the pore space of the intact rock and through the fractures is slow enough so that the convective heat transfer terms can be neglected. It is also assumed that the deformations of the intact rock and the fractures along with the fluid flow processes do not result in generation of heat. Therefore, heat transfer in the entire rock mass is described by the Fourier law of heat conduction: where k c is the equivalent thermal conductivity of the rock and h i are the components of the heat flux vector. Stress equilibrium equations written in terms of displacements u i are given by where n is the porosity, and ρ s is the density of the solid phase. From the fluid mass conservation law and Darcy's law (Eq. 3) the fluid flow equation can be derived as where α f is the thermal expansion coefficient of the fluid phase and K f is the bulk modulus of the fluid phase. From the thermal energy balance equation and Fourier's law (Eq. 4) the heat conduction equation can be obtained as where c p is the equivalent specific heat of the rock. Equations (5)-(7) constitute a system of governing equations for a coupled THM problem. If the effect of temperature is neglected, i.e. T ≡ 0, we recover the equations applicable to the coupled HM problem. Each of the Eqs. (5)-(7) is valid inside a domain with spatially uniform material properties, i.e. within the fractures or the intact part of the rock. The coupled THM problem described by the governing Eqs. (5)- (7) is  7) is solved using the Heat Transfer Module. By default these Modules are uncoupled, but the user can connect them by adding additional terms to the appropriate equations of each module as suggested by Eqs. (5) and (6).

The initial boundary value problem
Glacial loads on the surface of a rock mass arise due to the accumulation of snow, ice and water (Aschwanden and Blatter, 2009). The water content within a glacier is generally small (approximately 5 g of water per 1 kg of mixture) and the highest water content is found in the vicinity of the interface with the bedrock. The pattern of fluid flow both within and at the base of a glacier can be a complex THM process by itself. Problems associated with the mechanics of glaciers (Weertman, 1972;Marshall, 2005), the motion of glaciers (Fowler, 1981(Fowler, , 2011Morland, 1987;Picasso et al., 2004;Gomez et al., 2013;Stucki and Schlunegger, 2013;Ahlkrona et al., 2013;Thoma et al., 2014), the constitutive properties of ice and ice-rock interfaces (Picasso et al., 2004;Marshall, 2005;Headley et al., 2012;Tezaur et al., 2015), the contact with the bedrock and the evolution of their shapes (dell'Isola and Hutter, 1998;Picasso et al., 2004;Dietrich et al., 2010;Headley et al., 2012;Nielsen et al., 2012;Pollard and de-Conto, 2012;Gomez et al., 2013;Stucki and Schlunegger, 2013;Steffen et al., 2006;de Boer et al., 2014) have been studied. The glacier dynamics on a centennial timescale due to climatic warming was studied by Hambrey et al. (2005). Consider a model of a rock mass in the form of a parallelepiped with upper surface z = z top and lower surface z = z bot . Assume that a glacier of a non-uniform thickness moves along the upper surface of the rock, z = z top , parallel to the (x, y) plane. The thickness of the glacier is denoted by H = H (x, y, t), where (x, y) are the in-plane coordinates fixed at the base of the rock mass, as indicated in Fig. 1. As an approximation of the complex shape of the glacier, the profile of the glacier is often considered to be parabolic Steffen et al., 2014). In this case, the maximum ice-sheet thickness is expected to occur at the glacier centre and zero thickness at the glacier front. At the same time, it should be noted that the parabolic profile is very steep near the glacier front and very flat at the centre.
In order to examine the response of a poroelastic rock mass subjected to glacial loads, it is important to adequately take into account the interaction of the glacier with the underlying bedrock. On the temporal scale of interest to the glacial loading of a rock mass containing fracture zones, a simpler representation of the glacier can be adopted. Firstly, mechanical loading due to the weight of the ice sheet needs to be prescribed (Bangtsson and Lund, 2008;Read, 2008). The mechanical loading corresponds to the changing ice-sheet thickness H (x, y, t) and is applied as a time-and space-varying normal stress at the upper surface of the rock mass Chan and Stanchell, 2008;Read, 2008;Steffen et al., 2006;Person et al., 2012;Pollard and deConto, 2012). Since the densities of the glacier constituents, water and ice, are almost equal, the total normal stress acting on the upper surface due to the presence of the glacier is given by Secondly, hydraulic boundary conditions must be prescribed. Boulton et al. (1995), Boulton and Caban (1995),  and Person et al. (2012) state that if the permeability of the underlying rock is sufficiently low to discharge all basal meltwater as groundwater alone, groundwater heads are approximately equal to the ice pressures and the effective pressures are low. For example, Person et al. (2012) prescribed head boundary conditions on the upper surface of the sedimentary basin assuming that the head is equal to 90 % of the ice sheet elevation. In this study the intact rock is assumed to have a very low permeability and thus the fluid pressure acting on the upper surface can be approximated by Therefore, the effective stress exerted by the glacier on the upper surface of the rock mass can be evaluated from Eqs. (2) and (9) as To account for thermal loading caused by the presence of the glacier, the temperature change can be prescribed on the upper surface as where T u is a prescribed function.
On the lower surface of the rock mass, z = z bot , we prescribe zero values for the vertical displacement u 3 , the fluid velocity and the heat flux: On the vertical perimeter edges of the rock model we impose zero displacement, fluid velocity, and heat flux in the direction of the normal to the surface: where ∂/∂n = n∇ is the derivative in the direction normal to a surface. Initial fields within the rock mass correspond to the fields existing at the onset of glacier motion and thus should include stresses due to the weight of the glacier at its initial position, geostatic stresses, hydrostatic fluid pressure, the temperature distribution due to the geothermal heat flux, etc. In this work we focus our attention only on the change in the response of the rock mass caused by a small glacier advance and thus zero initial conditions can be prescribed, if the problem is linear. Figure 1 shows geometry of the fractured rock mass model used in the present study. The region of the rock mass is a parallelepiped with in-plane dimensions 22 000 m×22 000 m and a thickness of 1700 m. The fracture network is located at the interior of the rock mass domain, remote from its boundaries. The specific fracture network examined here contains a system of 21 planar fractures within a region of (approximately) 11 550 m × 10 200 m. Details of the fracture network are shown in Figs. 1, 2 and 10. Each dominant fracture has the shape of a prism or prismatoid. The thickness of each fracture is set to 10 m. This is a realistic size for a fracture; an example of a 10 m wide fracture zone situated in the Rhone Valley of the Alps, the Simplon normal fault, is given in the paper by Stucki and Schlunegger (2013). In the fracture network, 14 fractures are subvertical with a high dip angle (close to 90 • ), and 7 fractures are subhorizontal with a low dip angle (close to 0 • ). Among the 14 subvertical fractures, 13 fractures span the entire thickness of the rock model. Fractures that do not span the entire thickness of the rock (e.g. subhorizontal) originate at the upper surface of the rock. The subhorizontal fractures are numbered from 1 to 7 in Fig. 10.

Model geometry
Each fracture zone is defined by its surface which constitutes a planar quadrilateral formed by the four vertices lying in the same plane. To facilitate construction of the fracture surface, a 2-D workplane is constructed for each fracture. Orientation of the workplane is defined by specifying coordinates of the vector normal to the fracture surface. The fracture contour is drawn in the workplane by specifying coordinates of the four vertices of the fracture surface, transformed to the local coordinate system of the workplane. Subsequently, the fracture surface is extruded at a distance equal to the thickness of the fracture zone (10 m), forming a 3-D prism. Then this prism is embedded into a specified place within the geometry of the rock.
If a fracture originates at the upper surface of the rock, the upper edge of the fracture zone is expected to coincide with the upper surface of the rock, which is horizontal. However, after the extrusion process described above, if the fracture zone is not strictly vertical, the upper edge of the fracture (being orthogonal to the fracture surface) will be inclined with respect to the horizontal surface of the rock. Such a configuration may cause problems later during mesh generation since the edge of the fracture lies too close to the surface of the rock. To correct this geometry, the fractures' surfaces in their respective workplanes can be extended in such a way that, after extrusion and embedding, the fracture intersects the upper surface of the rock and extends beyond it. Then we simply cut off the portion of the fracture zone above the upper surface of the rock by using, for example, the subtraction operation. Similar corrections can be performed with respect to other edges of fractures that lie too close to the lower surface of the rock or to the perimeter surfaces. It should be noted that the fracture obtained in this way is a prismatoid and no longer an extruded feature.
Typically, it is desirable to have the capability of inserting into the 3-D geometry as many fractures as possible. On the other hand, if the geometry becomes too complex, the mesh may not be generated. Problems with mesh generation may arise when (a) the fractures are not parallel, i.e. have various orientations in space, (b) the number of intersections between fractures is too large, (c) the distance between the neighbouring fractures is too small, i.e. they are too close to each other, and (d) the thickness of fractures is too small. These conditions naturally place a limit on the number and the type of fractures that can be successfully inserted. We succeeded in inserting only 21 fracture zones into the geometry of the given rock mass as shown in Fig. 1. The basal rigid stratum will have an influence on the THM processes that are investigated. The depth of 1.7 km is chosen to conform to previous studies conducted by Chan and Stanchell (2008) so that there is a basis for comparison. When appropriate computing resources are available, the extent of the domain can of course be increased to accommodate depths similar to those employed by Ustazewski et al. (2008) when examining composite faults in the Swiss Alps formed by the interplay of tectonics, gravitation and postglacial rebound occurring at depth regimes of up to 5 km.

Mesh generation
The COMSOL Multiphysics® finite element software has a powerful mesh generator engine. For 3-D analysis only 2172 A. P. S. Selvadurai et al.: Thermo-hydro-mechanical processes in fractured rock formations tetrahedral elements can be used and by default the order of the elements is quadratic. There are several predefined element sizes such as "Extremely Coarse", "Extra Coarse", "Coarser", "Coarse", "Normal", "Fine", "Finer", "Extra Fine", "Extremely Fine" and, in addition, the user can define the element size, "Custom" mesh size. Each mesh is characterized by the following parameters: maximum element size, minimum element size, maximum element growth, resolution of curvature and resolution of narrow regions. The values of these parameters are geometrydependent. For example, for the given geometry of the rock mass, a parallelepiped with in-plane dimensions of 22 km × 22 km and a thickness of 1.7 km, the "Normal" mesh parameters have the following values: maximum element size 2290 m, minimum element size 412 m, maximum element growth 1.5, resolution of curvature 0.6, resolution of narrow regions 0.5.
If the model geometry is too complex, it might be impossible to construct the mesh with the predefined mesh size due to the error messages. For example, for the present geometry of the fractured rock mass, COMSOL Multiphysics® failed to construct the "Normal" mesh, but managed to construct the "Extra Fine" mesh consisting of 3 321 534 elements and a "Fine" mesh with 223 530 elements. To reduce the number of elements even further we define the "Custom" mesh size for which we assign the following values to the parameters: maximum element size 1690 m, minimum element size 172 m, maximum element growth 1.6, resolution of curvature 0.7, resolution of narrow regions 0.5. The resulting mesh has 158 728 tetrahedral elements and is shown in Fig. 2.
If the error message appears during mesh generation, it is recommended that the mesh size parameters be modified; even a small change in a parameter value can sometimes allow COMSOL Multiphysics® to successfully construct the mesh. For example, for the given geometry with 21 fractures, one could construct the "Custom" mesh if the minimum element size is set to 152-154, 157, 159, 161, 167, 168, 170, 172, 186 or 189 m in the range of all values from 152 to 192 m. Otherwise, a change in the geometry should be attempted, such as deletion of the fracture that caused the error or a small translation of this fracture in the x-y plane. For example, deletion of the sub-horizontal fracture 2, marked in Fig. 10, extends the list of minimum element sizes for which the "Custom" mesh can be constructed to: 152-154, 157, 160, 161, 163, 164, 168, 170, 172, 173, 175, 177, 178, 180, 189 or 190 in the same range from 152 to 192 m. We note that even if a certain fracture was removed from the model due to the error messages, another fracture can still be inserted.

Glacial loading
The primary emphasis of this work is to examine the THM response of the fractured rock mass when it is subjected to advancing glacial loads. In general, advance or retreat of the glacier can be caused by (a) sliding of the glacier along the bedrock and (b) a change in the accumulation and/or ablation rates. In the first case the shape of the glacier does not change, while in the latter case the shape of the glacier, i.e. its length and height, evolves. In the following we consider the glacier advance caused only by sliding, i.e. case (a). Figure 1 shows the glacier profile at time t = 0 and at time t > 0, after the glacier has moved a distance of 7500 m parallel to the x axis. The glacier front is assumed to be parallel to the y axis. Initially, at time t = 0, the glacier partially covers the model domain over the length of 7500 m; its front is located at x = −3159 m (Figs. 1, 3). It should be noted that the total glacier length is much larger than depicted in Fig. 1 but here only the marginal zone of the glacier that falls within the relatively small (compared to glacier length) dimensions of the model domain is shown.
The thickness of the glacier is non-uniform and can be described by the function H (X, t), where X is the distance from the centre of the glacier (datum point). Let H max the maximum thickness of the glacier (at the centre), and L be the glacier length, i.e. distance from the glacier centre to its front. Then the function H (X) at time t = 0 satisfies H (0) = H max and H (L) = 0. The specific form of the function H (X, t) adopted from the paper by Paterson (1972) is shown in Fig. 3 and described in Appendix A. The maximum glacier thickness is set equal to H max = 3200 m and H 2 max /L = 7.7, which gives the length L = 1 329 870 m. The profile of the ice sheet, shown in Fig. 3, is similar to that shown in Fig. 6 of the paper by Boulton and Caban (1995), in which, at a distance of 4000 m from the ice sheet front, the glacier thickness becomes equal to 300 m. The x coordinate of the model domain, shown in Figs. 1 and 3, and the glacier X coordinate are related by x = X − 1333029 m. The model domain boundaries are shown in Fig. 3 with vertical dashed lines.
The velocity of the glacier V is assumed to be 4.0306 × 10 −8 m s −1 , which is approximately 1.271 m year −1 . Thus, in 5900 years the ice sheet front will advance a distance of 7500 m across the model domain (see Fig. 3). This value of glacier velocity is about 100 times smaller than the one used by  and about 10 times smaller than the velocity calculated by Aschwanden and Blatter (2009). Such a small glacier velocity is chosen to avoid problems with a loss of accuracy of the solution, which arise when the glacier velocity becomes higher and the mesh is not sufficiently fine (see Sect. 4.1). In addition to the actual glacier profile, shown with a solid line, Fig. 3 shows an approximation to this profile with a dashed line. This approximation coincides with the actual glacier profile away from the ice sheet front; however, close to the leading edge of the ice sheet it replaces the actual profile with a much smoother function. In actual numerical computations, the use of this approximation is advantageous since it does not have any sharp corners; smoothing the loading profile has an extremely favourable effect on the convergence properties of the computational scheme, especially when using iterative solvers.
As mentioned, the initial or in situ state is not modelled in the present study and therefore the loading applied to the surface of the rock mass only corresponds to the difference between the current position of the glacier, at time t, and its initial position at t = 0. (This is an admissible simplification of the problem when the problem is linear.) In this case, the maximum glacier loading will always be applied at the initial position of the glacier front; the amount of the glacier loading will decrease away from the initial positioning of the leading edge.
In the following we focus our attention on the response of the rock mass over a short time (e.g. 5900 years), for a small glacier advance, when the glacier front is still close to the fractured region of the rock and the glacier centre remains at a very large distance from the fractured region. Thus, we do not compute the results corresponding to the maximum value of the glacier loads, when the glacier centre itself overruns the model domain. The choice of the time frame can be explained by the highly non-uniform distribution of the glacier load near the front of the glacier; our intention is to study the effect of this non-uniform load. In addition, we believe that the glacial sliding is limited to small distances compared to the overall glacier length.

Finite element solver
The response of the rock mass was obtained by using the transient finite element (FE) solver in COMSOL Multi-physics®. The numerical integration in time was implicit by default, which implies that a matrix factorization is required at each time step. Selecting the proper FE solver is of outmost importance because of the size of memory required, computational time limitations and issues related to convergence. In order to solve the HM problem (Eqs. 5-6) for the "Custom" mesh, a direct solver, such as the SPOOLES (SParse Object Oriented Linear Equations Solver) solver, requires a very large memory (larger than 30 GB) and a significant amount of time for matrix factorization. Thus, a direct solver cannot be used for the present problem. After testing several iterative solvers available in COMSOL Multiphysics®, it was found that the iterative GMRES (Generalized Minimal RESidual method) solver with the geometric multigrid preconditioner was, in fact, the only solver capable of dealing with this problem within a reasonable time frame (approximately 3 h of runtime when solving the HM problem and 8 h for the THM problem). On the memory side, GMRES required only slightly more than 6 GB of RAM. Other preconditioners had slower runtimes due to convergence issues.

Material properties
The fracture zones have THM properties substantially different from those of the intact rock. In particular, the permeability of the fracture zones is taken to be about six orders of magnitude greater than the permeability of the intact rock, whereas Young's modulus is assumed to be approximately 10 times smaller.
The properties of the intact rock and the fractures used in our computational simulations are the following: Young's modulus E D is 30 GPa for the intact rock and 3 GPa for the fractures; Poisson's ratio ν D is 0.25; the Biot coefficient α is 0.7; permeability k is 1.55 × 10 −19 m 2 for the intact rock and 3 × 10 −13 m 2 for the fractures; dynamic viscosity µ is 0.001 Pa s; fluid bulk modulus K f is 2.2 GPa; thermal expansion coefficient for the solid phase α s is 8.3 × 10 −6 1/ • C and for the fluid phase α f is 69 × 10 −6 1/ • C; equivalent specific heat c p is 1.8305 × 10 6 N/(m 2 K); and equivalent conductivity k c is 3.66 W/(m K) both for the intact rock and for the fracture zones. The porosity is assumed to be 0.05 for both the fracture zones and the intact rock. These material properties for the intact rock and fracture zones are similar to those 2174 A. P. S. Selvadurai et al.: Thermo-hydro-mechanical processes in fractured rock formations used by Chan and Stanchell (2008) to model the Canadian Shield subregion. Alterations to fracture permeabilities during glacial loadings are important but these extensions are not considered in this study. The estimation of the permeability of fractures in particular is governed by the hydraulic aperture of the fracture, which can be influenced by the normal and shear stresses acting on the fracture and gouge development (Nguyen and Selvadurai, 1998;Selvadurai and Nguyen, 1999;Selvadurai and Yu, 2005;Selvadurai, 2015). Discussions of these topics and references to further studies of this important topic are given by Boulon et al. (1993), Nguyen and Selvadurai (1998) and Selvadurai (2015).

Model performance
In this section we assess our choice of the user-defined "Custom" mesh, consisting of 158 728 tetrahedral elements (Fig. 2), and show that this mesh is suitable for solving the given problem of the bedrock response if the glacier speed is as low as 1.27 m year −1 . We will refer to this "Custom" mesh as "fracture-biased" mesh.

Mesh refinement evaluation for homogeneous rock
The first evaluation procedure involves analysis of the HM solutions for the homogeneous rock, i.e. the rock in which properties of the fractures are the same as those of the intact part. Here we examine how the presence of the densely refined regions in the mesh, in close proximity to the "fracture" zones, may affect the solution. To do so, in addition to the "Custom" mesh we consider the regular mesh with only 5786 elements that does not account for the presence of "fracture" zones. Figure 4 shows the distribution of vertical fluid velocity v z in the homogeneous rock for the two meshes described above. The time is 5900 years, during which time the glacier has moved a distance of 7500 m. The current position of the glacier front corresponds to transitioning from negative to positive velocity values, which is clearly seen in the figure. The magnitude of the negative velocity is considerably larger than the positive velocity. It can be seen that the fluid velocity distributions for the two meshes are quite similar, i.e. the presence of densely-meshed fracture zones does not change the solution significantly. This is an expected result since the rock is homogeneous. Our numerical tests also showed that the solution for the "Custom" fracture-biased mesh deviates significantly from the regular solution once the glacier speed is increased to 12.7 m year −1 . This implies that for a larger glacier velocity, the finer mesh will be required.

Mesh refinement evaluation for the hydraulic problem
The second evaluation procedure employed involves a comparison of the solutions for two meshes: the user-defined "Custom" mesh with 158 728 elements and the pre-defined "Extra Fine" mesh with 3 321 534 elements. Taking into account the fact that the "Extra fine" mesh has a massive number of elements, which would require a considerable amount of memory, we could solve only the "pure" hydraulic problem (Eq. 6) in which the deformations of the porous skeleton are suppressed, thus reducing the number of unknowns. Figure 5 shows the distribution of vertical fluid velocity v z in the intact rock for the two meshes described above. It can be seen that the vertical fluid velocity distributions for the two meshes are very similar. As before, the fluid velocity behind the glacier front is negative and it is positive ahead of it. The fluid velocity in the intact rock between the fractures is significantly reduced as indicated by the dark blue colour of the fracture-free region of the rock changing to green, yellow and even red in the region between the fractures. ] in the intact part of the fractured rock for "Custom" mesh (top) and "Extra Fine" mesh (bottom). Time is 5900 years and glacier velocity is 1.27 m year −1 . The velocity is obtained by solving a pure hydraulic problem.

Solution of HM problem for "Custom" mesh
Figures 6-15 show the hydro-mechanical response of the fractured rock subjected to glacial loading, described previously, at time 5900 years. The time-dependent fluid pressure and vertical compressive stress induced by the glacial load are applied to the upper surface of the rock mass as prescribed by Eqs. (8) and (9). The solution is obtained by solving the hydro-mechanical problem (Eqs. 5, 6) on the userdefined "Custom" mesh. For more accurate results, the model domain, shown in Figs. 1 and 2, was augmented along the x axis by two blocks of 22 000 m length each behind the glacier front, and by one block of 14 000 m length ahead of the glacier front. The total length of the augmented domain is therefore 80 000 m, while the size along the y axis remains the same at 22 000 m. Mesh size parameters retain the same values as for the user-defined "Custom" mesh, described above. The new blocks are fracture-free and are discretized using regular, non-biased meshes. The resulting mesh for the extended domain has 167 303 tetrahedral elements, just slightly more than the fracture-biased mesh for the main block of the domain, containing fractures. The increased length of the model domain enables us to take into account the glacier loading more accurately. Figure 6 shows the typical fluid pressure distribution at the interior, fractured region of the rock and exterior, homogeneous region at time 5900 years. To capture the response in the fractured region, the fluid pressure was plotted in the cross section where the y coordinate was equal to 6200 m, which is approximately at the mid-section of the rock mass (see Fig. 1). As is expected, the fluid pressure close to the lower surface is generally smaller than the pressure applied at the upper surface. The maximum fluid pressure occurs close to the initial location of the ice sheet front, i.e. to x = −3159 m. Within and around the fractures, the fluid pressure is almost uniform through the thickness as the pressure induced at the lower surface is almost equal to the pressure applied at the upper surface. Such a response is caused by the high permeability of the fractures. However, away from the fracture zones, the pressure at the lower surface remains smaller than that at the upper surface, due to the low permeability of the intact rock. This pressure difference creates a vertical pressure gradient ∂p/∂z and a fluid flow in the downward direction. For the purpose of comparison, Fig. 6 also shows the distribution of fluid pressure in the homogeneous region of the rock in a cross section parallel to the x axis; the influence of the fractures in altering the fluid pressure distribution in the region due to glacier loading is quite evident. In particular, we can observe that the pressure in the intact rock in close proximity to the fractures is larger than the pressure in the fracture-free, or homogeneous, region of the rock, for the same value of the x coordinate. This implies that the vertical fluid velocity in the intact rock near the fractures is generally expected to be smaller than that in the homogeneous region of the rock (see also Figs. 7 and 8). Figure 7 shows the fluid pressure variation at a depth of 200 m from the upper surface plotted as a function of the x coordinate for two cross sections: y = 6200 m, in the fractured region of the rock mass, and y = −3200 m, in the homogeneous region. The time is 5900 years. As before, the fluid pressure and stress changes caused by the glacial advance of 7500 m are applied at the upper surface of the rock mass. From the figure it is evident that the fluid pressure in the intact part of the fractured rock is larger than that in the homogeneous region of the rock, for the same value of the x coordinate. This implies that the pressure gradient and, therefore, fluid velocity in the intact part of the fractured rock is smaller, in absolute value, than that in the homogeneous rock (given that the fluid pressure applied at the upper surface slightly exceeds 2 MPa).   The vertical velocity distribution in the rock is complex. The fluid velocity can be both positive and negative. A negative fluid velocity implies that fluid flow takes place into the rock (i.e. downwards) and a positive fluid velocity implies that the fluid flows out of the rock (i.e. upwards). It can be seen that, in general, the magnitude of the negative velocity is much larger than the positive velocity, which is consistent with the direction of the meltwater flow. In fact, beneath the glacier most of the fluid flows into the rock with a velocity of −3 × 10 −13 m s −1 and, exterior to the glacial footprint, most of the fluid flows out of the rock with a velocity 2.5 × 10 −14 m s −1 . Between the fractures the magnitude of the negative fluid velocity in the intact rock is in general smaller compared to that in the homogeneous region of the rock (as the dark blue colour in the exterior region of the rock changes to light blue, yellow and even red in the interior region at close proximity to the fractures). Figure 9 shows variation of the vertical fluid velocity in the intact part of the rock mass with fractures in the cross sections y = 6200 and −3200 m corresponding to the fractured and homogeneous regions, respectively, at a depth of 200 m from the upper surface, at time 5900 years. The figure clearly shows how the presence of fractures affects the fluid velocity in the intact rock near the fractures: the fluid velocity in the intact rock of the fractured region is smaller, in absolute value, than the velocity in the fracture-free region of the rock. Figure 10 illustrates the vertical fluid velocity v z within the fractures of the rock at 5900 years. It can be seen that in the fractures oriented subvertically (most fractures are subvertical) the fluid velocity is much larger than that in the intact part of the rock (compare Fig. 10 with Fig. 8). In fact, the maximum fluid velocity in the fractures can reach a value of −5.5×10 −8 m s −1 , which is about 50 000 times larger than in the intact rock, but still within the Darcy regime. The greater fluid velocity in the fractures can be explained by the large permeability of the fractures compared to that of the intact rock. We also observe that the fluid can flow both in the negative and positive directions even within a single fracture that is inclined to its strike with respect to the glacier front, i.e. not parallel to the glacier front. For example, the largest subvertical fracture in the model is inclined to the strike at an angle of 34.7 • with respect to the y axis. At one end of this fracture,  the vertical fluid velocity is towards the upper surface of the rock and can reach the value of +5 × 10 −8 m s −1 , although this fracture is totally covered by the glacier at 5900 years. The fact that some fluid in the fracture can flow towards the upper surface can be explained by the non-uniform distribution of the ice sheet loading along the x axis and the very high permeability of fractures. It should also be noted that in the fractures with a low dip angle (subhorizontal orientation) the vertical fluid velocity is generally smaller than that in the subvertical fractures. Figure 11 shows variation of the vertical fluid velocity in the largest fracture zone along its length, at a depth of 200 m  from the upper surface, at time 5900 years. At this depth the velocity distribution is approximately uniform across the thickness of the fracture. It is apparent from the figure that at this depth the velocity in the fracture remains high and negative at one end of the fracture (−3×10 −8 m s −1 ) and positive at another end (+5 × 10 −8 m s −1 ). Figure 12 shows the vertical displacement or deflection in the rock mass with fracture zones at time 5900 years. The displacement is negative under the glacier, which is indicative of the mechanical compression of the rock, and positive outside the glacier, which is indicative of the conventional heave associated with glacial loading (Walcott, 1970(Walcott, , 1976McNutt and Menard, 1978;Selvadurai, 1979Selvadurai, , 2009bSelvadurai, , 2012 Selvadurai and Nguyen, 1995). The maximum negative displacement (−0.076 m) is near the front of the glacier at its initial position. It can be observed that the displacements in the fractured and homogeneous regions of the rock are not very different, i.e. the displacement field is not strongly influenced by the presence of fractures with low deformability properties. Figure 13 shows variation of the vertical displacement in the rock mass with fracture zones in the cross sections of the fractured and fracture-free regions y = 6200 and −3200 m, respectively, at a depth of 200 m below the upper surface, at time 5900 years. The vertical displacement is in fact smaller, in absolute value, in the fractured region of the rock although the difference between the displacements in the fractured and fracture-free regions of the rock is not significant. Figure 14 shows the mean effective stress (σ xx + σ yy + σ zz )/3 in the intact part of the rock with fractures at 5900 years. The mean effective stress is mostly compressive but tensile stresses do occur in the fractured region at the upper surface (dark red patch). The maximum compressive stress occurs near the lower surface of the rock close to the initial location of the ice sheet front. The distinctive feature of the response of the rock with fractures is the existence of tensile stresses near the upper surface between and around the fractures (dark red patch); however, the magnitude of these tensile stresses is smaller than the compressive stress. It can also be observed that in the homogeneous, exterior region of the rock, the stress is entirely compressive. It should be noted that for any failure analysis the HM stress state resulting from glacial loading should be combined with the in situ geostatic stress, GIA stress (Steffen et al., 2006) and possible tectonic background stresses.   Figure 15 shows the mean effective stress variation in the intact part of the rock with fracture zones in the cross sections y = 6200 and −3200 m corresponding to the fractured and fracture-free regions of the rock, respectively, at a depth of 200 m below the upper surface, at time 5900 years. At this depth the mean effective stress in the fractured region of the rock remains tensile and can reach the value of +0.1 MPa whereas in the fracture-free region this stress is compressive, approximately equal to −0.2 MPa.

Solution of THM problem for "Custom" mesh
In the thermo-hydro-mechanical problem the temperature change caused by glacier advance or growth must be prescribed at the upper surface of the rock.  examined the change in the ground surface temperature through the entire glacial cycle using climate/temperature simulation. From this simulation it was clear that before the arrival of the glacier, the ground surface temperature decreases significantly from +10 to −25 • C. As the glacier gradually occupies the domain of interest, the insulating effect of the ice sheet gradually increases the subglacial temperature from −25 to 0 • C. This complex temperature evolution at the upper surface of the rock mass is simplified in the present work and only a net cool down from +10 to 0 • C, associated with the glaciation period, is considered. Similar surface temperature evolution was used in the study by Read (2008): the mean annual surface temperature was set to +7 • C during non-glacial periods and, to account for glaciation, the mean annual rock mass surface temperature was reduced to 0 • C.
In Figs. 16-19 we illustrate the response of the rock mass with fracture zones subjected to the glacial load described previously and, in addition, to the effects of cooling associated with the advance of the glacier. The time-dependent fluid pressure, the vertical compressive stress and the temperature change are applied to the upper surface of the rock mass as prescribed by Eqs. (8), (9) and (11). As indicated previously, the initial or in situ state is not taken into account in the present study and, therefore, only the surface loading and temperature that correspond to the difference between the current position of the glacier, at time t, and its initial position are applied to the surface of the rock. The temperature change on the upper surface caused by glacier advance is assumed to have the following form: where H (x, t) is the smoothed/approximated glacier thickness at time t (Appendix A). The maximum temperature change is set equal to T max = −10 • C and will be reached at the initial location of the ice sheet front when t = L/V = 1 047 000 years, i.e. when the centre of the glacier overruns the model domain. Since in this study attention is restricted only to a small glacier advance, this time can never be reached during computational simulations, i.e. the glacier motion will be terminated a long time before this time is reached. Figure 16 shows the temperature change generated within the rock 5900 years after initiation of the glacier motion. The temperature change on the upper surface is prescribed according to Eq. (14). The temperature change under the glacier reaches the maximum of −0.66 • C at the upper surface of the rock close to the initial location of the ice sheet front and is almost zero beyond the glacier margin. (The current maximum of −0.66 • C is much smaller than T max = Figure 16. Temperature change distribution [ • C] within the rock mass 5900 years after the start of the glacier advance. −10 • C since the centre of the glacier is very far from the model domain boundaries at time 5900 years.) We observe that the distribution of the temperature through the thickness of the rock mass under the glacier is non-uniform, which implies that the temperature did not reach the steady state in 5900 years and at the lower surface of the rock the temperature is still above the value prescribed at the upper surface. Figure 17 shows the variation of the temperature change at the upper surface of the rock mass and at a depth of 800 m from the upper surface, at time 5900 years. The dependence of the temperature field on the y coordinate is only due to a difference in the mesh refinement in the fractured and homogeneous regions and thus rather weak. From the figure it is evident that the temperature at a depth of 800 m is much smaller than the prescribed temperature at the upper surface, which suggests that even for such a small glacier speed a steady state was not reached in 5900 years. Figure 18 shows the mean effective stress in the intact part of the rock with fractures at time 5900 years. The rock mass is subjected to the compressive stress, the fluid pressure induced by the glacial load and the temperature field (Eq. 14) associated with the glacial advance. We note that the compressive mean stress at the lower surface of the rock mass is almost equal to that in the rock subjected only to the mechanical load and fluid pressure (see Fig. 14). However, the tensile mean stress, which occurs in the fractured region in the vicinity of the upper surface, is increased from 0.15 to 0.28 MPa. In contrast to the effective stress, the vertical fluid velocity is not significantly influenced by the cooling; the velocity field in the rock subjected to the compressive stress, the fluid pressure and the temperature field (Eq. 14) is almost the same as that in the rock subjected only to the mechanical load and fluid pressure (as shown in Figs. 8 and 9). Figure 19 shows variation of the mean effective stress in the intact part of the rock mass with fractures at a depth of 200 m below the upper surface in the cross sections  through the fractured and fracture-free regions y = 6200 and −3200 m, respectively, at time 5900 years. At this depth the mean effective stress in the rock of the fracture-free region is compressive, approximately equal to −0.15 MPa, whereas the stress in the intact rock of the fractured region remains tensile, reaching the value of +0.15 MPa.

Concluding remarks
In the present work, the behaviour of a fractured rock mass subjected to advancing glacial loads was analysed. The glacial load was taken into account by prescribing a spatially and time-dependent compressive stress and fluid pres- sure on the upper surface of the rock mass. Resulting fluid pressures, vertical fluid velocities, displacements, and mean stress distribution diagrams were obtained through computational simulations that incorporated THM processes. Results were plotted for the specific time of 5900 years during which the glacier advanced across the model domain by a distance of approximately 7500 m.
The computational results show that there is a difference between the fluid pressures at the upper and lower surfaces of the rock mass due to the low permeability of the intact rock. This pressure difference creates a vertical pressure gradient that contributes to the vertical fluid velocity component, which is a major feature of the fluid flow in the present problem. We also observed a difference in the fluid pressure distributions within the fractured region and the surrounding homogeneous region of the rock. In fact, the fluid pressure within the fracture zones and the intact rock in close proximity to the fractures is larger than that in the homogeneous region, located remotely from the fracture zones. This can be attributed to the large permeability of the fractures.
The distribution of the vertical fluid velocity within the intact part of the fractured rock is complex. In general, the direction of the fluid velocity is negative (downward) beneath the glacier and positive (upward) beyond the glacier margin. The magnitude of the downward velocity is larger than the upward velocity, which is consistent with the direction associated with the meltwater flow. As was already noted, between the subvertical fractures and near horizontal fractures the vertical fluid velocity in the intact rock is generally smaller than that in the homogeneous region. The vertical fluid velocity within the fractures of the rock is about 50 000 times larger than the velocity in the intact rock. The direction of the fluid velocity in the fractures is mainly negative but large velocities in the positive direction can also occur, especially in the subvertical fractures that are not parallel to the ice sheet front. The role of the fracture zones and their dominant influence on the fluid flow in the region during a glaciation event is thus confirmed by the present study.
The mean effective stress in the intact part of the rock mass is mostly compressive, but tensile stresses do occur between and around the fractures, at close proximity to the upper surface of the rock. The maximum compressive stress is at the lower surface of the rock at the initial position of the ice sheet front. The glacial loading induces compression of the rock, i.e. negative strain in the direction of the rock thickness. For the situations examined in the paper, this axial compression is relatively uninfluenced by the presence of fracture zones. In addition to the compressive stress and fluid pressure, the given rock mass was subjected to the cooling associated with a glacier advance. Since the present study considers a small glacier advance of only 7500 m, the associated temperature change is also small, only −0.66 • C. It was observed that even for this small temperature change, the mean effective tensile stress noticeably increases at the upper surface of the rock between and around the fractures. However, the effect of the cooling on the vertical fluid velocity distribution in the fractured rock was small. These results suggest that the hydro-thermal interaction in the present problem can be neglected, but the effect of cooling on the stress distribution in the intact part of the rock (thermo-mechanical interaction) should be given greater consideration.
This study is closely related to the work performed earlier by  and Stanchell (2005, 2008) in that it also examines the effect of the glacial loading on the behaviour of fractured bedrock. In their work the rock was modelled as a thermo-poroelastic medium and contained 3-D (but narrow) fracture regions, which were also assumed thermo-poroelastic. The novelty of the present model is the use of a quite complicated fracture network consisting of 21 fractures whereas only 5 fractures were included into the models described in the papers by  and . The NWMO report by Chan and Stanchell (2008) deals with a more complex fracture network consisting of 19 fractures but the fracture thickness is set to 20 m, whereas in the present work it is 10 m, which implies greater geometry and mesh complexity. In addition, in the present study the focus is on the transient response caused by highly non-uniform and moving glacial loading acting on the surface of the bedrock. This loading is directly linked to the specified glacier velocity and profile. Another novelty of this work is the investigation of the use of a commercially available finite element code, which can reduce the efforts of the users in solving the present problem by finite element method.
Computational modelling of fractured rock regions capable of accounting for THM processes is an essential tool for examining the response of rock masses at timescales relevant to glaciation scenarios. The computational approach allows for the examination of the relative importance of components in the THM model and the study of their influence on both the intact rock and the fracture zones at timescales relevant to glaciation scenarios and ultimately to deep geologic disposal concepts. In the modelling considered in this paper, both the fracture zones network and the intact rock are assumed to be sessile. The modelling can be extended to include activation of faults, e.g. in a combined analysis with large-scale models of GIA; natural hydraulic heterogeneities of the geologic media Selvadurai, 2010, 2014), particularly, at close proximity to the fracture surfaces; and stress-induced alterations of intact permeability. These effects can have an important influence on groundwater flow patterns in the domain. The information on the influences of stress states on permeability evolution in fractures is particularly scant and in a detailed analysis relevant to a deep geologic repository setting this aspect needs to be considered. Also implicit in the study is the absence of both radiogenic and geothermal heating effects; the influences of the latter could be important on the timescale of a glaciation event. The initial position of the glacier can be described by the function H (X), where H is the thickness of the glacier and X is the distance from the centre of the glacier (the datum point). The function H (X) satisfies H (0) = H max and H (L) = 0, where H max is maximum thickness of the glacier and L is the distance from the centre to the glacier front.
The glacier has accumulation and ablation zones. For simplicity, we assume that the net thickness of ice added to the surface each year in the accumulation zone is equal to the net thickness of the ice lost each year in the ablation zone. Owing to this assumption, the accumulation and ablation zones have equal lengths, L/2. In this case, the profile of the glacier in the steady state can be obtained as (Paterson, 1972) H (X) = 2 1/8 H max L − X L 1/2 , in the ablation zone L/2 < X < L; H (X) H max 8/3 + 2 1/3 X L 4/3 = 1, in the accumulation zone 0 < X < L/2.
Knowing the velocity of the ice sheet V and the initial profile, given by (A1), we can obtain the equation for the glacier profile in the ablation zone as a function of time: H (X, t) = 2 1/8 H max L − X + V t L 1/2 , in the ablation zone L/2 < X − V t < L.
The approximation to the parabolic profile (Eq. A2) is evaluated by the formula The source files used in the computations performed in connection with this paper can be found at the following links: https://www.dropbox.com/s/83b1gcrukw47sib/ glacierproblem_mesh158728_noheat.mph?dl=0, https: //www.dropbox.com/s/oy7fo28303td1lf/glacierproblem_ mesh158728_withheat.mph?dl=0.