DECOVALEX-THMC Project. TASK E. Implications of Glaciation and Coupled Thermohydromechanical Processes on Shield Flow System Evolution and Performance Assessment. Final Report

A Deep Geologic Repository (DGR) situated on the Canadian Shield will be subject to long-term climate change that will markedly alter surface conditions as a result of glaciation and permafrost penetration. Task E of DECOVALEX THMC is a systematic numerical case study of the subsurface THM processes and mechanisms arising from long-term climate change. The case study focuses on predicting the magnitudes and rates of change in groundwater flow and state of stress caused by time-varying glacial boundary conditions acting on a 1.6-km deep, subregional scale (∼100 km2), fractured Shield flow system adapted from Sykes et al. and includes a simplified version of Realization 1 of the stochastic Fracture-zone Network Model (FNM) of Srivastava. A suite of transient, 2D and 3D MOTIF finite-element, coupled THM, subsurface simulations were completed. Depth-dependent fluid salinity and temperature-dependent fluid density and viscosity have been included, as well as initial and transient thermal, hydraulic and mechanical boundary conditions developed from the 'nn2008 (base case) and nn2778 (alternative case) realizations' of the Univ. of Toronto's GSM model for two 120,000-year Laurentide glaciation scenarios. The GSM model reveals that for either scenario, the Shield subregion in question was subjected to three cycles of ice-sheet advance and retreat, with varying magnitudes and rates of change in basal normal stress, temperature, meltwater fluxes and permafrost evolution. Key findings of Task E coupled THM subsurface modelling include: i) for this particular conceptual model, the flow domain appears to have little memory of previous glacial cycles, with respect to carryover of significant thermal and hydraulic effects; ii) the increase of hydraulic head under ice loading, primarily caused by consolidation effects, is not equal to the total stress imposed by the glacier on the bedrock, but rather is about 1/3 of the glacial basal normal stress, in part due to the ratio between the compressibilities of the rock and the water; iii) the incremental head values resulting from ice loading and HM coupling are relatively uniform throughout the modelled subregion both horizontally and vertically, in contrast to what would occur if one used a hydraulic head boundary condition equivalent to the ice-sheet thickness in an uncoupled flow model; iv) head values in fracture zones (FZ's) differ from those in the adjacent rock mass (RM) by a few metres; v) during the short-lived glacial advance and retreat stages, the horizontal hydraulic gradient is steeper than the nonglacial topographically driven hydraulic gradient (∼0.003) by up to a factor of two, while near the Last Glacial Maximum (∼20 ka BP) the hydraulic gradient is slightly less than the nonglacial value; vi) throughout the glacial cycle the vertical head gradient is generally <0.05 from surface to ∼725m below surface, while below this level, saline groundwater with an assumed density of 1.10 kg/L gives rise to a vertical hydraulic gradient of ∼0.1; vii) sensitivity analyses showed that a combination of a temperate glacier, very low permeability rock (∼10-20 m2 ) and limited FZ connectivity is necessary for residual anomalous hydraulic heads to persist at depth for thousands of years following deglaciation; viii) effectively truncating most of the FZ interconnectivity by means of 2D modelling yielded small anomalous heads on the order of 10m; ix) including density effects from infinity was found to result in calculated anomalous heads, mainly hydrostatic heads in saline groundwater, being substantially greater at depth x) during the glacial cycle, Darcy fluxes (velocities) in the FZ's and highly permeable RM near surface are on the order of 10-2 to 10-1 m/a and range between 10-7 to 10-5 m/a in the RM below 350 m; xi) within the context of the Task E conceptual model, FZ sets with different orientations or belonging to different interconnectivity groupings exhibited relatively consistent groundwater velocities; xii) through using 2D simulations and thereby removing most of the FZ interconnectivity, groundwater velocities in the FZ's were reduced by a factor of 100; xiii) conservative particle-tracking analysis indicated that 72% of glacial meltwater particles did not penetrate more than 70m below surface and only 6% penetrated to 500m or further below surface; xiv) meltwater penetration depths are slightly enhanced by thermal effects, slightly diminished by density effects from depth dependent salinity, slightly enhanced by simulating a smooth glacial topography, hardly influenced by the glacial scenario (within the two scenarios simulated), but severely underestimated by using a 2D model that truncates most of the FZ connectivity; xv) during glacial advance/retreat, principal effective stresses are reoriented, factor of safety is slightly enhanced in the RM and reduced in the FZ's; xvi) salinity, glacial topography and the glacial scenario have little effect on the predicted effective state of stress, but 2D modelling drastically overestimates the stability of the FZ's. A limited 2D numerical study was conducted on subregional Shield groundwater flow dynamics under permafrost conditions. With a time-invariant, uniform 350-m thick permafrost layer modelled as low-permeability (∼10-19 m2 ) rock, coupled HM simulations with salinity predicted low Darcy flux in the permafrost, but somewhat higher flux in the RM unit immediately below, compared to when no permafrost was present