The reactivation of inherited tectonic structures formed during the Paleoproterozoic Trans-Hudson Orogeny (THO) has played a significant role in generating high-grade unconformity-related uranium deposits in the eastern Athabasca Basin. The role of these tectonic structures is now investigated through a series of two-dimensional hydrothermal numerical models. Two modelling scenarios are considered: (1) models during the THO peak of metamorphism and (2) models with a permeable layer mimicking the presence of the Athabasca Basin, deposited unconformably over the THO basement. In the first scenario, general fluid patterns are strongly affected by the applied permeability configurations. Unidirectional high fluid flow zones (from 10-9 to 10-8 m·s-1) and high thermal gradients (up to 65°C·km-1) can be observed above and within the deep-seated tectonic structures. In the second scenario, well-established fluid convection cells or unidirectional fluid flow zones are observed within the basin layer, with upflow originating from the core of the deep-seated structures, regardless of the applied fluid pressure regime. These results highlight that these deep-seated structures can efficiently transport fluids and heat towards the upper parts of the crust and the basin. In the second scenario, the loci for preore alteration are then evaluated by computing a rock alteration index based on temperature and fluid velocity constraints. These alteration areas reside along and above the deep-seated structures and are potential regions for structural reactivation during mineralization. These results imply that the analysis of the inherited tectonic structures, combined with the alteration regions, can serve as markers for uranium exploration.

1. Introduction

The formation of mineral deposits requires multiple processes related to hydrothermal systems. According to the “mineral system” framework, four critical elements must be defined to evaluate the conditions to generate mineral deposits: lithosphere architecture, favourable geodynamics and transient self-organized hydrothermal systems, the fertility of the geological domains, and the preservation of the depositional environment [1, 2]. This framework indicates that large long-lived, deeply penetrating, and steeply dipping structures that juxtapose distinctly different basement domains are geological features that lead towards giant deposits [1]. Several studies have evaluated how the “mineral system” framework can be applied from the regional to the camp scale and translated in terms of project generation and exploration targeting [1, 3, 4]. This concept is currently guiding geophysical surveys, highlighting the limits of crustal and lithospheric blocks [5, 6], the reassessments of the potential of regional metallogenic provinces (e.g., [3, 7]), and the detailed mineralogical and geochemical studies of the successive hydrothermal alteration episodes [812]. However, the contribution of numerical modelling to constrain the role of such deeply penetrating and steeply dipping structures as long-lived vectors for fluids and heat transport has not been investigated.

The Athabasca Basin (Saskatchewan, Canada), a Paleo- to Mesoproterozoic siliciclastic basin, hosts giant unconformity-related uranium deposits and a major metallogenic province [13]. At the eastern extent of the Athabasca Basin, multiple high-grade and high tonnage unconformity-related uranium deposits reside along a NE-SW trending basement-hosted structure that extend over several hundred kilometres along strike. This is the Wollaston-Mudjatik Transition Zone (WMTZ) formed during the Trans-Hudson Orogeny (THO; Figure 1(a); [14, 15]). As revealed by exploration drill holes and its magnetic geophysical signature below the Athabasca Basin, this tectonic structure is composed of steeply dipping anastomosed shear zones (Figure 1(b), [1618]). At its western extent, recent discoveries of basement hosted world class uranium deposits along the Patterson Lake conductive trend [19] and extensive geophysical and geological surveys [6, 2022] reveal the presence of crustal to lithospheric scale deeply rooted structures. These structural features can be considered as critical elements for the location of major hydrothermal systems and generation of giant uranium deposits.

At the unconformity-related uranium deposits in the eastern Athabasca Basin, uranium mineralization has been related to the steeply dipping inherited basement structures, considering their ability to move large volumes of fluid (e.g., [2326]). However, the mechanism for fluid transport and metal deposition remains unclear. Li et al. [25] argued that short-term deformation-driven flows might be incapable of transporting sufficient metalliferous-enriched fluid needed to form large unconformity-related uranium deposits. They proposed that thermally driven fluid convection is more capable of transporting large volumes of uranium-bearing fluids to deposition sites and that alternating deformation-driven fluid flow and thermally driven fluid convection are responsible for uranium mineralization. The establishment of convection cells at thermal equilibrium is long-lasting, even though fault reactivations play an important role in enhancing the permeability of the conduits.

This numerical modelling study is aimed at addressing the following questions concerning fluid flow patterns and physical conditions required to form uranium deposits within and above deep-seated structures during periods of quiescence: (1)What were the fluid flow patterns and their main drivers during the metamorphic retrogression of the Wollaston-Mudjatik Transition Zone?(2)What was the role of the ancient deep-seated inherited structures and how the pre-Athabasca fluid flow pattern and related alteration impacted or controlled the location of the hydrothermal activity at the origin of the unconformity-related uranium deposits?(3)What are the implications of deep-seated structures for future ore exploration?

For this study, we consider two 2-dimensional modelling scenarios that include inherited tectonic structures formed before the formation of the Athabasca Basin. The first modelling scenario mimics the tectonic conditions following the evolution of the WMTZ during the peak of metamorphism of the THO. The performed models allow to test different permeability configurations of increasing complexity before the deposition of the Athabasca Basin. The second modelling scenario involves a permeable layer representing the presence of the Athabasca Basin. A rock alteration index constrained by fluid velocity and temperature distribution [27] is then computed to predict zones where preore chlorite alteration is most likely to occur when the second modelling scenario has reached steady-state conditions.

The crustal architecture of the eastern part of the Athabasca Basin is inherited from the Paleoproterozoic Trans-Hudson Orogeny (THO; Figure 1) which resulted from the collision between the Hearne craton, the Sask craton, and the Superior craton between ca 1.87 and 1.79 Ga [15, 2834]. The Hearne craton, consisting of Archean granulitic tonalitic-trondhjemitic gneiss domes of the Mudjatik Domain and the Paleoproterozoic metasedimentary gneisses of the Wollaston Domain, is interpreted as a passive margin sequence [31]. These domains were intensively imbricated and deformed during the Trans-Hudson Orogen that led to the formation of the WMTZ [14, 28, 32]. The WMTZ is characterized as subvertical anastomosing fault systems that are highlighted by narrow magnetic low signatures (Figures 1 and 2(a); [1618, 35]). Mylonite zones are characterized by steeply dipping foliations bearing subvertical stretching lineation [36, 37] and by the presence of Paleoproterozoic sediments that reached amphibolite metamorphic facies conditions (i.e.,  GPa, °C, [14, 15, 38]).

Previous thermomechanical numerical models have verified that such deep-seated major shear zones could be formed from a protracted shortening of the lithosphere [39] and that sediments can reach burial - conditions in agreement with the peak of metamorphism known for the THO at ca. 1.84–1.82 Ga ([39] and Figure 2(b)). Thermobarometry studies indicate a decrease in pressure from 0.8–1 GPa to ~0.5 GPa occurring at 1.82 to 1.77 Ga due to exhumation of the Archean basement and changes in the tectonic regime [15, 38]. Subsequently, a protracted cooling period occurred in the region at 1.76–1.71 Ga, marking the late-THO tectonic evolution [28, 38, 4042].

Tectonic evolution of the THO is finally sealed by the deposition of detrital sediments that form the Paleoproterozoic to Mesoproterozoic basins (including Athabasca) which rest unconformably over the exhumed domains, preserving the underlying basement-hosted structures. These steeply dipping tectonic structures within the basement and in particular the graphite-bearing shear zones were reactivated after the deposition of the Athabasca Basin. This reactivation is marked in many fault segments by an offset of the unconformity that reaches about 100 m in the McArthur River uranium deposits [13, 16] and by the development of cataclastic breccia and gouge in the basement faults and brittle faults into the basin (e.g., [13, 4347]).

This polyphase tectonic evolution and the existence of favourable transient events are also considered as critical elements to characterize the unconformity-related uranium mineral system [16]. The continuity between the reactivated faults of the basement and the brittle faults in the basin and the development of alteration haloes and mineralization within these structures suggest that the reactivation has controlled the circulation of the mineralizing fluids ([13, 16, 48]). Geochemical and fluid inclusion studies indicate that uranium mineralization is coeval with the tectonic reactivation of deep-seated inherited structures [26, 4952]. Previous hydrothermal numerical results indicate that mineralization can be associated with tectonic reactivation occurring after the thermal blanketing effect generated by the Athabasca Basin [24, 25, 5355]. The way how fluids circulate between the basement and basin through the inherited structures leads to two types of unconformity-related uranium deposits characterized by different alteration zones [13, 49, 5658]. The first type is the sandstone-hosted type associated with the discharge (egress) of the fluids from the basement. The sandstone-hosted type deposits have two end-member alteration patterns that are located within the sandstone overlying the unconformity: (1) zones of quartz dissolution+illite and quartz overgrowth+illite-kaolinite-chlorite+dravite as observed at the Cigar Lake deposit [43] and (2) zones of silicified-dravite-chlorite-kaolinite and illite-dravite as characterized in the McArthur River deposit [46]. The second type is the basement-hosted deposits resulting from recharge (ingress) of the fluids from the basin. Narrow, inverted alteration halos of sudoite±illite and chlorite+biotite±sudoite are located at the outer margins of the inherited basement structures [58]. However, the spatial distribution of the alteration halos is extensive, up to 400 m wide at the base of the sandstone, and it may exceed several thousand metres in strike length (e.g., [13]). These unconformity-related uranium deposits and their alteration zones attest to the widespread fluid-rock interactions and fluid circulation between the major fault systems of the basement and the Athabasca Basin.

3. Modelling Approach

3.1. Governing Equations

The hydrothermal models are performed with FLAC3D (the Fast Lagrangian Analysis of Continua in 3 Dimensions). This code solves the coupled fluid flow and heat transfer equations on a Cartesian grid via an iterative process that involves sequential stepping of fluid and thermal modules [59]. It has been widely used for modelling fluid flow, heat transfer, and deformation, especially at an outcrop to upper crustal scale [24, 5355, 6063].

Heat transfer is solved by the following series of equations that define the energy balance in a convective-diffusive heat setting. Conservation of energy (Equation (1)) is achieved through where is temperature; is the thermal flux; is the fluid specific discharge; is the volumetric heat source intensity; and are the reference density and specific heat capacity of the fluid, respectively; and represents the specific heat capacity of the rocks.

Conductive heat transfer (Equation (2)) is solved by Fourier’s law which defines specific heat flow as heat flow normalised by area: where is the thermal conductivity of rocks and is the temperature gradient vector.

The temperature values from the previous time step (Equation (3)) are used to calculate the fluid density in the following time steps, which is expressed as where is the volumetric thermal expansion of the fluid at . The calculated fluid density is used to compute fluid velocity according to Darcy’s Law (Equation (4)), which is expressed as where is the fluid dynamic viscosity, is the gravitational acceleration, is the assigned permeability of the rock, and is the fluid gradient pressure.

Boussinesq approximation is applied to the models.

3.2. Model Setup
3.2.1. Geological Setting for Modelling Sections

The first main geological scenario, M1 (Figure 2(b)), mimics the evolution of the WMTZ during the peak of metamorphism of the THO. The model consists of two homogenous layers of rocks simulating the presence of metasediments (representing the Wollaston Supergroup) deposited unconformably on top of the Archean granitic basement. The geometry of the metasediments that intersects the basement is determined by the sedimentary cusps formed due to pop-down tectonics [36, 39, 64, 65]. From M1, three modelling subscenarios will be evaluated as M1A, M1B, and M1C (see Section 3.2.3 for more details).

The second main geological scenario, M2 (Figure 2(c)), takes into account a period of tectonic quiescence of 120–160 Ma occurring between the peak of metamorphism of the THO (∼1.84 Ga; [14]) and the deposition of the Athabasca Basin (~1.76–1.72 Ga; [38, 42]). During this period, the eroded region undergoes a slow postorogenic cooling [41]. The region underwent 20 km of distributed exhumation leading to an uplift of the basement up to the depth of 7 km. This erosion phase removed most of the complex structural features, leaving behind a remnant of the deep-seated structures starting at the depth of 6 km. The resultant geological model agrees with previous seismic surveys and thermotectonic interpretations [14, 66, 67]. From M2, three modelling subscenarios will be evaluated as M2A, M2B, and M2C (see Section 3.2.3 for more details).

3.2.2. Geometry, Boundary, and Gridding of the Models

The upper and lower boundaries are left open for fluids to leave the system, whereas the side boundaries remain closed and impermeable. Initial fluid pressure at the bottom boundary is set to the calculated fluid pressure (see Section 3.2.3 for more details). Initial temperature at the top of the model is maintained at 20°C. At the bottom boundary, the initial temperature is assigned with a value based on a linear geothermal gradient of 30°C·km-1.

All M1 models have a physical model domain of 200 km length, 0.1 km width, and 15 km depth (see Figures 3(a)–3(c)). For M1A and M1B in Figures 3(a) and 3(b), the structural model was designed using an equidistant rectangular grid with a cell resolution of 1 km by 0.15 km (20000 elements; Figure 4(a), A1 and A2). The width of the vertical structures penetrating the basement is 3 km. For M1C, a Tetrahedral Irregular Network (TIN) was used to generate a more accurate representation of the geometry for the structures and deformed regions (Figures 3(b) and 4(b)). With 17710 elements, high-resolution zones are assigned to where we expected to have high hydrothermal activity (i.e., within the metasediments and tectonic structures) to reduce computation time (Figure 4, B1 and B2). The size of the grid for low-permeability regions such as the basement is assigned a coarser grid (up to 2 km spacing in Figure 4, B2).

To focus our study on the hydrothermal activity involving the basin and the surviving deep-seated structures, the length of M2 models is reduced from 200 km to 110 km to provide focus in our study. The same TIN grid system of 14514 elements is also applied (Figure 4(c)). Between the depths of 3 to 8 km, the high-permeability zones (i.e., basin and deep-seated structures) are refined to have high-resolution regions (tetrahedral grid smaller than 200 m, seen in Figure 4, C1 and C2). Low-resolution grids (i.e., tetrahedral grid larger than 1.5 km) are employed to the basement and deformed basement rocks as both represent permeability regions.

3.2.3. Model Fluid Configuration

The hydrothermal modelling protocol follows a simple physical procedure that does not consider vapour-fluid multiphase, fluid-fluid, and fluid-rock interaction. Among the fluid properties, only density varies with temperature. Other parameters such as fluid viscosity, heat capacity, and thermal conductivity remain constant (see Table 1). According to Sibson [68], fluid pressure increases with depth following different regimes (see the blue dashed line in Figure 5). The upper part of the profile (from 0 to 4 km depth) reflects the hydrostatic pressure conditions where the pore fluid factor, (where and are the fluid pressure and the vertical stress, respectively). Below, the fluid pressure profile corresponds to the transition from a hydrostatic to a lithostatic regime. This regime (where ) can extend up to 15 km depth (see [68]). However, below 10 km depth, fluid pressure is very close to lithostatic pressure () without considering other potential horizontal stresses.

In our models, a simplified dual-pressure regime is considered (black dashed line in Figure 5), as seen in previous published numerical models (e.g., [69]). Hydrostatic pressure is first applied up to 7 km depth () in order to account for the porous sandstone-dominated nature of the Athabasca strata that form a relatively permeable environment [49, 50, 70, 71]. Fluid pressure increases to near lithostatic pressure from 7 to 8 km depth allowing to simulate the transition from hydrostatic to lithostatic pressure.

For the M2 models in Figure 4(c), three variants of fluid pressure regimes are considered. M2A uses the same simplified dual-pressure regime as in M1. The second model, M2B, follows the same simplified dual-pressure system of M2A but fluid pressure evaluated at the bottom of the model is kept constant. It allows to simulate the presence of a constant fluid flow maintained throughout the more permeable regions (i.e., deformed basement and deep-seated structures). M2C uses a simplified hydrostatic pressure regime only (i.e., ) to evaluate the effect of fluid pressure change in fluid flow patterns.

3.2.4. Model Permeability Configurations

Estimating rock permeability () in any geologic model is required for fluid flow modelling but assigning permeability values to rocks is challenging, considering their wide range of possible values and their sensitivity to physical processes (e.g., [72, 73]). Hence, a suitable permeability value is required primarily for large and deep-seated structures (10 km by 15 km) to establish a baseline of the fluid flow patterns. Previously published numerical works simulated large and deep fault zones with a constant permeability of 10-15 m2 that extends to depths of 15 km (e.g., [27, 61, 74]). High fluid flux values in our initial numerical tests of 10-15 m2 for the deep-seated structures cause fluid flow values not to converge, resulting in the model being numerically unstable. The permeability of the deep-seated structures (  m2) is then decreased to allow the fluid flow and temperature equations to converge, thus achieving a thermally steady-state condition.

In M1 subscenarios, permeability of the rocks is first varied to evaluate their effect on the fluid flow patterns. The porosity and permeability assigned to the various rocks remain constant regardless of model scenario (Table 1). Three permeability configurations are considered. The first subscenario (M1A in Figure 3(a)) has a bimodal permeability configuration with a of 10-16 m2 for the metasediments and a of 10−18 m2 for the basement. The second subscenario (M1B in Figure 3(b)) utilises a simplified depth-dependent permeability function for metasediments, as indicated by experimental data [7577]. The depth-dependent function is expressed as , where (m2) and (km) correspond to rock isotropic permeability and depth, respectively. The permeability of the basement is kept constant. High-strain regions found within the metasediments and basement are incorporated into the third subscenario (M1C in Figure 3(c)) to simulate tectonic inheritance (M1 in Figure 2(b)). The spacing and geometry of the inherited zones are then in good agreement with the anastomosed fault systems observed nature (in Figure 2(a)). These high-strain regions in the models for both metasediments and basement are labelled as deformed metasediments and deformed basement, respectively. The permeability of these deformed metasediments ( m2) and basement (  m2) is assumed to be more permeable than their nondeformed counterparts that mimic the damage zones during fault generation (e.g., [78, 79]). The change of permeability from their nondeformed counterparts is highly uncertain due to the heterogeneity of the fault zone and sampling process [79]. The permeabilities of the deformed metasediments and basement are considered constant.

The Athabasca Basin in M2 is composed of a 2 km permeable homogenous sandstone layer (  m2) followed by a 4 km low-permeability eroded cover ( m2) that represents the clay-rich impermeable Wolverine Point Formation overlying the Manitou Falls Formation (Figure 3(d); [42]). The thermal conductivity of the remaining deep-seated structures (in particular, the exhumed deformed metasediments cover located within the basement) was increased (from 2.5 to 5.0 W·K−1·m−1) to reflect the presence of alteration minerals (i.e., graphite, quartz, and chlorite) that have been identified as vectors for uranium uptake and precipitation [49, 55, 8083].

3.3. Thermally Constrained Restricted Rock Alteration Index (RRAI)

Hydrothermal fluids contain dissolved chemical species that flow through rock columns. Tectonic events that lead to strain localisation perturb the environment via permeability enhancement, which allows the transport of reactants in solution towards the reaction site by advection and diffusion (e.g., [73, 8486]). Potential alterations and mineralization patterns can be deduced by looking at zones where the cooling rate is essential. This cooling rate is referred to as the “rate of mineralisation” [87] that corresponds to the negative values of the rock alteration index (RAI). RAI is evaluated as , where and represent fluid velocity and temperature, respectively. RAI has been used in previous hydrothermal models to evaluate and predict ore deposition potential (e.g., [63, 74, 88]). RAI values highlight the importance of the deep-seated structures concerning unconformity-related uranium deposits. Negative RAI indicate fluid flows to regions of low temperature. Potential zones for mineralization and alteration can be predicted in a metamorphic setting when  m·s-1 (see [27], for more details). This fluid velocity value is then applied as a constraint to RAI which we now consider as the restricted rock alteration index (RRAI).

Alteration halos delineate the presence of these unconformity-related deposits following the main deposit types. Preore chlorite alteration is then considered as a pathfinder and a potential chemical trap for uranium mineralization [9, 28, 89]. We propose that RRAI can be constrained further by associating it with the formation temperature of the preore chlorite alteration. These temperatures are found to vary widely in deposits: from 275 to 300°C in the Millennium deposit [9]; 230–250°C at Dawn Lake, McArthur River, and Rabbit Lake [89]; and 235–245°C in McArthur River [11]. Therefore, in combination with RRAI and the formation temperature of preore alteration, this RRAI certifies that the physical conditions for chlorite alteration are met; if fluid temperature is within 230–300°C and RAI is negative and  m·s-1. This thermally constrained RRAI will only be applied to the M2 models because formation temperatures of preore alterations are found in deposits located during the presence of the Athabasca Basin and the metasediments of the WMTZ, and the earliest uranium mineralization event occurred with the presence of the Athabasca Basin (ca. 1.68 Ga; [28]).

3.4. Model Limitations

The primary assumption in our numerical models is that the fluid viscosity is kept constant. In FLAC3D, rock permeability interacts with the rock’s isotropic permeability coefficient and is evaluated as , where and are hydraulic conductivity and gravitational acceleration, respectively. Furthermore, intrinsic permeability () is evaluated to be inversely proportional to dynamic viscosity (). Making fluid viscosity variable will cause unwanted permeability enhancement to the initial permeability setup seen in Sheldon [61]. In the context of the unconformity-related uranium deposits at the Athabasca Basin, the source of the basement-associated fluids is highly debatable (e.g., [12]). Considering the complex geological layout from the thermomechanical model, we can only investigate the fluid flow patterns by keeping fluid viscosity constant and preventing any permeability enhancement. Variations in the fluid viscosity can be considered in future model iterations.

4. Results

4.1. Model M1

In the bimodal permeability configuration for M1A, macroscopic fluid flow patterns indicate unidirectional fluid flow within the metasediments. Fluid is observed to flow along the contact between the impermeable basement and in regions where sedimentary cusping occurred at  km, -20 km, 0 km, and ca. +25 km. No convection cells are observed (Figure 6(a)). Fluid temperature in the metasediments reached at least 500°C where sedimentary cusps are present and 475°C when they are not present at depth. Three regions of increased fluid flow velocity of at least an order of magnitude (-11 m·s-1 to 10-10 m·s-1) in comparison with the basement are located: (1) within the metasediments where sedimentary cusps are present at depth, (2) below the contact of the metasediments and the basement, and (3) within the basement at ~13 km depth. This region of increased fluid velocity corresponds to the transition from hydrostatic to lithostatic pressure (as indicated in Figure 5) and is present in all M1 models (Figure 6).

The presence of a variable permeability in metasediments of model M1B enhances a more variable and periodic fluid flow. Upward flows occur where sedimentary cusps are present (at  km, -20 km, 0 km, and +25 km in Figure 6(b)), and downward flows occur aside them (at  km, -30 km, +10 km, and +60 km in Figure 6(b)). The pairing of upward and downward flows coupled with the corresponding deflections of the isotherms indicates the presence of well-established convection cells.

Unidirectional fluid flow vectors in M1C (where in M1C is reduced by an order of magnitude) indicate that no convection cells form in these conditions (Figure 6(c)). This, compared with model M1A (Figure 6(a)), prevents the generation of high temperatures at the upper parts of the model. Fluid velocities show regions of high fluid velocities ( m·s-1) that correspond to the geometries of tectonic inheritance present in the basement and metasediments. Upward deflections of the isotherms indeed clearly define areas of increased fluid flow and geometries of tectonic inheritance. Additionally, fluid vectors are also observed to be oriented preferentially along these regions. At the transition between hydrostatic to lithostatic pressure regime, the temperature is at least 250°C (ca. 40°C km-1) and up to 400°C (ca. 70°C km-1) at the undeformed basement and the opening of the deep-seated structures, respectively.

4.2. Model M2

In model M2A, the primary convective cells occur directly above the deep-seated structures, with a wavelength of 5 km and high fluid velocities of 10-7 m·s-1. Temperatures at the contact between deep-seated structures and the basin peaked at 350°C and 375°C (see Figure 7, A1 and A2, respectively). At the 6 km depth, a temperature of 375°C in the basin is near the contact of the deep-seated structure. Between the 6 and 8 km depths, regions of increased fluid velocity ( m·s-1 in Figure 7(a)) are generated at the contact of the deformed basement and deep-seated structures and the impermeable basement. This depth range also corresponds to the depth of the transition from hydrostatic to lithostatic fluid pressure (in Figure 5), which results to higher fluid pressure within more permeable regions. Fluid flows within the basin layer indicate the presence of well-established convection cells (in Figures 7(a) and 7(b)). Widespread convection cells within the basin follow downward deflections of the 200°C isotherms. The wavelengths of these convection cells within the basin are smaller than the convection cells above the deep-seated structures (from 2.5 km to 1 km). Fluid temperature in these convection cells is between 200 and 225°C. Around the basin, fluid velocity is around 10-11 to 10-10 m·s-1. When compared to M2A, fluid flow patterns, velocities, and temperature distribution in M2B are similar even though the applied dual fluid pressure regime in M2B has an open boundary (in Figure 7(b), B1 and B2). It is apparent that fluid velocity is controlled by the permeability of the rocks as the permeability profile is similar in the two models.

Regions of increased fluid pressure (i.e., 6–8 km and 12–14 km depth) are not present in model M2C (Figure 7(c)) unlike in M2A and M2B. Well-established convection cells are present at the contact of the deep-seated structures, with a wavelength of 5 km and fluid velocities of 10-7 m·s-1. At the contact between the deep-seated structures and the basin, the elevated temperature peaked at 300°C at 6 km depth. Fluid flows near the deep-seated structures indicate two types of flow within the basin layer (in Figure 7, C1 and C2): (1) a lateral flow at the lower parts of the basin layer and (2) a vertical fluid flow, as seen as seen from the cyclical deflection of the 200°C isotherm at the upper part of the basin layer (Figure 7, C1 and C2). The wavelength and amplitude of the 200°C isotherm deflections are small (<500 m). Furthermore, fluid flow in the basin close to the deformed basement ( to +10 km in Figure 7(c)) presents some well-established convective cells with short wavelengths (<500 m).

Loci for chlorite alteration can be found in three areas in M2A and M2B (in Figures 8(a) and 8(b)). Alteration region 1 is located at ca. 7 km depth and corresponds to high fluid velocity area within the basement. It marks the transition from hydrostatic to lithostatic pressures (in Figure 5). This alteration region corresponds to basement-hosted and not confined by the deep-seated structures. Alteration region 2 corresponds to zones associated with the deformed basement, bounded by the 300°C isotherms, and is located at 6 to 8 km depth. Here, fluids have the potential to flow from zones of lower to higher permeability [24]. The vertical geometry of the chlorite alteration aligns with the geometry of the deformed basement. The loci of the chlorite alteration is comparable with observations from Millennium deposit (e.g., Figure 3 in [9]). Alteration region 3 corresponds to the loci of chlorite alteration bounded by the 300°C isotherms within the basin layer and above the contact of the deep-seated structures. This locus of alteration region 3 is comparable to the distribution of the illite-dominated alteration region observed in Cigar Lake [43]. These results show that, at the first order, a dual fluid pressure regime can set the physical conditions allowing to create the alteration distribution for a basin-hosted unconformity-related uranium deposit type.

In M2C, the loci for chlorite alteration are restricted to alteration regions 2 and 3 (in Figure 8(c)). As it is bounded by the 200 and 300°C isotherms, alteration distribution within the deep-seated structures is limited. The loci of alteration in M2C in the basin are more extensive than M2A and M2B. The presence of well-established convective cells in our model may have sufficient fluid-rock interaction to form large alteration haloes as illite in Cigar Lake (e.g., [24, 50]).

5. Discussion

5.1. General Fluid Patterns and Convection Cells along and above Deep-Seated Structures

In all M1 models, fluid convection cells occurred only in M1B (Figure 6(b)) in a depth-dependent permeability setting. These convection cells occur due to the presence of more impermeable zones in the metasediments at depth and the cooling effect of the top boundary (20°C). Applying less permeable zones in higher temperature regions towards the bottom boundary reduces fluid velocity, allowing thermal conduction to dominate. The resultant temperature field becomes more homogenous when a less permeable metasediments zone in M1C is applied ( m2 in Figure 6(c)). Reducing the rate of fluid-thermal advection also prevents the cooling effect of the top boundary from getting overwhelmed, causing cooler and denser fluids to flow downwards. In M1A, the absence of the downward fluid flow indicates that the cooling effect from the top boundary is insufficient and is overwhelmed by the high fluid advection.

Two main drivers of fluid flow are present in all M2 models regardless of the applied pressure regime (Figure 7): (1) thermal conduction is dominant within an impermeable basement as fluid flow is less active. This leads to unidirectional fluid flow and horizontal isotherms. (2) Thermal convective flow is dominant in permeable rock regions (i.e., the deep-seated structures and the basin). This type of flow is characterized by the deflections of the isotherms, high fluid velocities, and the convective motion of fluids. These convection cells occur in all M2 models. The maximum temperature of the fluids flowing above the deep-seated structures in M2A and M2B (ca. 350°C in Figures 7(a) and 7(b)) is higher compared to M2C (ca. 290°C in Figure 7(c)). This increase in temperature of the fluids is caused by the presence of a fluid pressure gradient from lithostatic (at depth) to hydrostatic (towards the surface) within the deep-seated structures. This fluid pressure gradient leads to an increased fluid flux and, consequently, to a heating controlled by advection.

In all M2 models (in Figure 7), widespread convective cells within the basin layer are consistent with previously published numerical results (e.g., [24, 25, 54]). The main difference is the scale of the numerical model. Li et al. [25, 54] considered fluid flow activity on a local deposit (i.e., the McArthur River deposit at the 5 km by 7 km depth), whereas this study considers fluid flow activity at a regional scale (i.e., WMTZ at 150 km by 15 km depth). Additionally, the deep and large structures in the M2 models are considered as large homogenous rock units within an impermeable basement (i.e., 3–5 km wide). In contrast, the numerical models conducted by Eldursi et al. [24] and Li et al. [25] considered the effects of one and multiple discrete structures mimicking superficial fault systems at a deposit scale in a homogenous layout, respectively. Our results in the M2 models show that well-established thermal convections do occur above the deep-seated structures. At the deposit scale, our results indicate preferential fluid flow within these structures in the basin and in the deeper parts of the basement (see Figure 4 in [24]). The depth of this fluid flow within the structures also corresponds to the transition from hydrostatic to lithostatic pressure regime, favouring advective flow [90]. Our results indicate that the presence of these deep-seated structures is essential to initiate thermal convections. Furthermore, both occurrences of long-lived convection cells above the deep-seated structures and fluid flow occurring preferentially within these structures also explain why the Wollaston-Mudjatik Transition zone is more favourable for mineral deposits.

5.2. Role of Inherited Deep-Seated Structures in Preparing the Geological Environment for Mineralization and Significance for Exploration

The depth where the transition from lithostatic to hydrostatic pressure regime occurs fulfils one of the conditions for the permeable structures to operate as potential fluid-pressure-activated valves [68]. These fault systems are indeed favourable loci to form alteration zones and correspond to preferential area where local deformation can occur during a transient tectonic events [14, 29]. Furthermore, the reactivation of these deep-seated structures can cause localised fluid-thermal fluxes, creating microscale networks of permeability swarms and deposit minerals (e.g., [8, 26, 91]). At a local deposit scale, hydrothermal fluids do flow out and into the structures during compressional and extensional strain, respectively [24, 25, 53]. The onset of fault activation can also temporary disrupt prevailing thermal convection within the basin [25].

Preore chlorite alteration, which acts as a chemical reductant in proximity to graphite, has been identified as a potential chemical trap ([28] and references therein) and as an essential vector for uranium mineralization in the Millennium deposit [9] and P-Patch [52] and Kiggavik-Andrew Lake structural trend prospects (e.g., [92]). The deep-seated structures are efficient conduits for transporting fluids and heat from the basement to the basin. Thus, they become the ideal loci for hosting key reductant minerals such as graphite and rehydrate metamorphic minerals such as chlorite (e.g., [13, 52]). Hydrothermal simulation results indicate that the fault systems maintain some localised fluid-thermal fluxes (i.e., natural hydrothermal state). They also show that these structures have the physical conditions for chlorite alteration to occur prior tectonic reactivation regardless of permeability configuration and pressure profile (see Figures 7 and 8, respectively).

Dynamic pressure on these structures may also play a vital role in fault reactivation. Hence, total pressure within the brittle crust is estimated to be up to ~2 times the lithostatic pressure under compression and ~2/3 times the lithostatic pressure under extension [9396]. These brittle structures would become focal points for pressure variation, depending on their rheology (e.g., [97, 98]). In the Athabasca Basin, it is postulated that the impermeable seal of the fault and fluids is subjected to this pressure variation due to far-field stresses. Sources for these far-field stresses could be from the deposition of the thick (~7 km) Athabasca Basin itself (e.g., [42]), protracted deep erosion and exhumation of the THO (e.g., [41]), and the formation and breakup of supercontinents Columbia/Nuna and Rodinia (e.g., [99]). The oldest preore alteration and uranium mineralization in the Athabasca Basin are ca. 1.68 Ga and ca. 1.59 Ga, respectively [28] which provides about 90 Myrs for tectonic loading to occur. At the onset of fault formation, permeability increases significantly and is followed by high fluid flow from depth to the basin layer. The amount of overpressured fluids ejected from fault failure can be significant, especially during the first few years (e.g., [100]). Hydrothermal precipitation occurs if the physical and chemical conditions for preore chlorite alteration and uranium mineralization are met (Figure 8). This is also in good agreement with the fault-valve hypothesis proposed by Sibson [86, 90]. However, this appears to be a transient process. After all, hydrothermal fluids have the potential to precipitate additional minerals such as quartz, uraninite, and chlorite within the newly made permeability swarm and change the physical conditions needed for rock failure [8, 91, 96].

At a regional scale, magnetotelluric profiles have been conducted to characterize the structure of the lithosphere and structures’ association to mineral deposits (e.g., [5, 6, 101]). Magnetotelluric profiles conducted towards southern Trans-Hudson Orogeny by Jones et al. [5] and in the Taltson orogeny towards the SW Athabasca Basin show low-resistivity anomalies that have a close spatial relationship with high-strain tectonic domains and to major mineral occurrences [6, 21, 22]. These observed anomalies extend over strike for hundreds of kilometres, (e.g., North American Central Plans, [5]). Moreover, field data, such as the presence of carbonatites within the associated deep-seated shear zones of the Patterson Lake corridor, indicates that these deep-seated structures reach mantle depths [20]. Similar crustal-scale low-resistivity corridor networks are also well documented in the Gawler craton, Australia, and refer these high-strain zones for potential fluid flow [22, 102104]. Subsequently, these tectonic structures are assumed to be strongly modified through long-standing chemical and physical processes. The rocks associated with these deep-seated inherited structures in our M2 models would have undergone multiple syntectonic depositions, tectonic reactivations, and fluid flow activity for at least 120–160 Myrs (i.e., from the peak of metamorphism of the THO to the deposition of the Athabasca Basin). Fluid flows preferentially from the impermeable basement to the permeable basin through the deep-seated structures to provide fluids and heat, potentially fulfilling critical conditions for uranium mineralization [13, 23]. Such interpretations correspond very well to the tectonic structures generated from the thermomechanical models proposed in Poh et al. [39, 65] and the structural cross-sections ([17, 35] and Figure 2(a)).

Therefore, the role of the deep-seated structures can be surmised in two main points: (1) natural candidates for fluid flow channelling and fluid-rock interaction and (2) preferential zones for tectonic reactivation during transient geodynamic activity.

5.3. Refinement of Exploration Strategies

The unconformity-related uranium deposits hosted by the Athabasca Basin are spatially associated with three main tectonic zones that still have the potential for further exploration: (1) WMTZ, (2) Snowbird tectonic zone, and (3) Patterson Lake deformation corridor. These deep-seated structures can establish convection cells with a significant temperature increase (up to 350°C) and discharge a significant amount of fluid from the basement to the basin (Figure 7). Physical conditions within the basin, deep-seated structures, and the deformed basement are indeed suitable for potential preore chlorite alteration in our models and are in good agreement with previous fluid inclusion and clay mineral studies ([9, 57, 82, 83, 105, 106]; Figure 8). The predicted alteration regions also represent a high potential for a chemical trap located within the basin layer, the deep-seated structures, and the deformed basin (Figure 8). Naturally, these structures become potential sites for fluid accumulation and high fluid-rock interaction. Within the WMTZ, the footprint of unconformity-related uranium ore bodies is spatially associated with second-order fault zones or transecting late brittle faults ([16, 17], and Figure 2(a)). It is suggested that the exploration strategy can be expanded to include deep-seated structures at depth and out of the basin, and it will potentially lead to significant discoveries.

6. Conclusions

This work is carried out on two tectonic scenarios mimicking the WMTZ before and after the Athabasca Basin formation and using 2-D hydrothermal models. The aim is to establish the role of deep-seated structures on fluid flow and the consequences to mineralization potential. Results show that deep-seated structures (up to 30 km depth) constitute fluid-thermal conduits capable of bringing fluids and heat towards the upper portions of the crust. Additionally, well-established and robust convective cells can occur near the deep-seated structures. Loci for preore chlorite alteration are also determined when applying the RRAI. The distribution of these areas strongly depend on the geometry and permeability of the inherited deep-seated structures. Furthermore, our results also show that these preore chlorite alteration corridors within the deep-seated structures also indicate that mechanically weak zones are prone to fault reactivation. This study may provide insight into the potential hydrodynamic activity in regions presenting a basin cover. An upscaling from camp to regional scale to properly characterize these deep-rooted, ancient structures and paleo-plate boundaries is, therefore, a viable exploration strategy.

Data Availability

Data can be made available upon reasonable request.


This work represents the third and final part of the author’s Ph.D. thesis.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This study was supported financially by Orano Canada Inc (Project number 39M.PSGEPA0425). Numerical simulations were performed at the GeoFluids lab at the University of Regina. The authors would like to thank Thibault Duretz, Laurent Guillou-Frottier, Remy Chemillac, and the team of Orano for their feedback and discussions on earlier versions of the manuscript.