Abstract
Gold deposits are not uniformly distributed along major faults due to complex (and long-debated) interactions between seismicity, hydrothermalism, and structural heterogeneities. Here, we use static stress modelling (SSM) to quantitatively investigate these interactions, by exploring the role of Cadillac-Larder Lake Fault (CLLF) Archean seismicity in the genesis of the regional goldfields. Various rheological factors are evaluated for optimizing the models’ ability to reproduce known gold occurrences, regarded as the fossil primary markers of synkinematic hydrothermal systems. We propose that the marked structural heterogeneities of the CLLF induced persistent seismic segmentation and recurrent ruptures of the same fault windows that arrested on robust node points. These ruptures favour repeated occurrences of seismically triggered hydrothermalism along long-existing fluid pathways having an enhanced permeability and iterative ore formation into supracrustal discharge zones by means of episodic drops and build-ups of pressure. Two-dimensional SSM permits the predictive mapping of these high-potential zones. These modelled zones correlate positively with the actual observed gold distribution. We demonstrate that the ruptures along the Joannes Segment arresting on the Davidson Fault and Lapa’s bend can explain the occurrence and location of the Rouyn and Malartic goldfields; the models’ validity is improved by implementing regional geological constraints; and the distant gold occurrences from the CLLF, including the Bourlamaque field, can be explained by doublet seismic events along the Rivière-Héva and Lapause subsidiary faults. Our results provide new perspectives from a fundamental standpoint and for exploration purposes.
1. Introduction
Mesothermal gold deposits are developed by the discrete circulation of mineralized fluid flows that are controlled by structural frameworks (e.g., [1–5]). The channeling and concentrating of the fluids are important for gold endowment. The role of fault seismicity in such processes has been discussed over the last three decades [6–10]. Yet, quantifying the regional-scale interplay between seismicity, crustal permeability, and hydrothermal ore genesis remains an open question requiring multidisciplinary efforts that involve both active and fossil systems. Addressing these issues from a mechanical modelling perspective is particularly challenging in fossil systems due to the lack of quantifiable evidence. The complex distribution of gold deposits along regional-scale faults has long been a major focus for understanding the controls on mesothermal gold mineralization processes [11–13]. Recent studies have shown that the nonuniform distributions of gold endowments are correlated to variations in fault geometry, including deviations from its main trend [14–16]. Most faults can be divided into several segments characterized by constant attitude. The segments are separated from each other by bending step-overs, hard links, or gaps [17]. These large fault segments play a critical role in seismic behaviour, displacement, and failure propagation. Historically, segmentation has been a major focus of seismology for understanding seismic hazards around existing active fault zones. This has produced the aftershocks theory; this idea involves the geomechanical modelling of static stress changes around a ruptured segment for predictive mapping of coseismic activity and aftershock triggering zones. The zones are generally referred to as coseismic damage zones (San Andreas Fault, California [18–20]; North Anatolian Fault, Turkey [21]). The approach has been validated using modern active faults, and the modelled zones correspond to actual aftershock locations 65% to 85% of the time [20, 22, 23].
From a hydrothermal standpoint, damage zones undergo a temporary and dramatic increase in structural permeability related to the opening of fracture networks where aftershocks are triggered [6, 24]. This results in the activation of second-order faults that lead to the draining of deep-seated fluid reservoirs [25] and a greatly increased hydrothermal circulation [8, 10]. The repeated failure of the same fault segment (hundreds to thousands of times), combined with the cyclic draining of pressurized fluids, allows for the iterative enrichment of gold deposits by the repeated circulation of large volumes of metal-rich fluids within a restrained area [7, 15]. From the perspective of gold metallogeny and prospection, predicting the spatial distribution of these potentially auriferous coseismic discharge zones—by means of a static stress modelling approach as analogous to seismological studies on present-day active faults—is a highly valuable and relevant approach. Micklethwaite and Cox [15, 26] pioneered this approach by applying it to the analysis of gold deposition near auriferous Archaean faults. They established that the activation of favourably oriented segments along the Black Flag and Boulder-Lefroy faults, Western Australia, yielded coseismic damage zones that were positively correlated with gold distribution; these zones include the world-class St. Yves and Mt. Pleasant deposits. Moreover, Micklethwaite and Cox [26] suggested that static stress changes from the mainshock combined with doublet aftershocks enhance the precision when delineating damage zones. These studies have focused mainly on the impact of some intrinsic characteristics of the main fault such as its geometry [14], permeability [27], and links between segments [28]. However, it is likely that the regional geological background (structural framework and lithologies) in the vicinity of the main fault greatly interferes with the coseismic propagation of stress and strain. This, along with the influence of preexisting, low-strength regional fabrics on coseismic stress patterns, remains an unaddressed first-order parameter. Although they play an important role in the spatial distribution of damage zones, the trends of the primary failing planes at any point near the ruptured main segment have also been given little attention [20]. Moreover, the validation of static stress models using ancient fault systems represents a major challenge [29] due to the absence of direct aftershock records. On the other hand, mesothermal gold mineralization located near fault areas represents reliable markers of seismically triggered fossil hydrothermalism occurring within damage zones. This study thus follows the premise used in previous studies (e.g., [9, 15, 26]) that, through a classic inverse approach, the validity of various models can be evaluated by comparing the simulated coseismic damage zones to the spatial distribution of actual goldfields.
This study aims to identify the role of Archean seismicity along the Cadillac-Larder Lake Fault (CLLF) on gold deposition in the Abitibi greenstone belt, Québec, Canada. Numerical geomechanical modelling is used to establish static stress changes subsequent to mainshocks along the fault. Emphasis is placed on testing the incorporation of regional geological inputs as mechanical and rheological constraints for refining the models.
The CLLF is a 250 km long east-striking structure associated with several clusters of gold deposits, which represent 25% of all gold produced in Canada [5]. The Abitibi Subprovince is characterized by numerous anastomosing faults dominated by east- and southeast trends [35, 36]. Rabeau [37] performed pioneer geomechanical modelling attempts in the Rouyn-Noranda region, using a strain propagation approach in a 3D elastic space. Dextral-reverse finite-displacement was applied on a 25 km long portion of the CLLF and subsidiary Wasa Fault, combined to compressive far-field boundaries conditions. Several methodological aspects deeply differ from the present study: the input kinematics combine strain fields related to different tectonic events, with likely different faults deformational regimes, and fluid-focusing zones were tracked by mapping high bulk volumetric changes while the present study computes the stress fields applied on fractures in order to track structural permeability enhancement in damage zones. Rafini [38] modelled damage zones of the CLLF during the late dextral strike-slip. The present study builds upon this latter modelling attempt. A century of gold exploration in the region has yielded considerable knowledge regarding the geological framework, including the architecture of major structures and lithologies. Therefore, this environment is well suited for testing the impact of contrasting lithologies and subsidiary faults on the geometry of stress change following displacement along a main fault such as the CLLF.
2. Geological Setting and Timing Relationships between Deformation and Gold Mineralization
2.1. Regional Geology
The Abitibi Subprovince is an Archean greenstone belt dominated by volcanosedimentary rocks and sedimentary basins intruded by several types of granitoid rocks (Figure 1). The Southern Volcanic Zone of the Abitibi Subprovince [39] is separated from the Pontiac Subprovince to the south by the CLLF. The Pontiac Subprovince is mainly composed of sedimentary and granitoid rocks. The CLLF can be separated into several 30–40 km long segments of constant attitude despite having geometric variations (azimuth, dip) between them (see [40, 41] for details). Tips between segments include abrupt and gentle bends as well as hard links such as second-order faults (Figure 1).
(a)
(b)
The Abitibi and Pontiac Subprovinces have experienced prolonged deformational events and are characterized by highly developed structural features. Rock assemblages are tilted prior to ductile deformation events [40, 42]. Although there is no consensus on the interpretation of the structural evolution of the study area, it is commonly accepted that a major N–S shortening event was responsible for most of the finite deformational framework. This was followed by a dextral strike-slip along east-striking faults during a subsequent NW–SE shortening (e.g., [40, 41, 43–45]). Most of the contacts between rock groups have been activated as first- and second-order faults [35, 42], as illustrated in Figure 2. Rocks are affected by a predominantly east-striking vertical foliation related to shortening. Regional metamorphism is mostly at a greenschist grade, but ranges locally from subgreenschist to amphibolite facies [46, 47]. This metamorphism grade is consistent with the lower levels of a seismogenic crust [6]. The exposed portion of the CLLF experienced a frictional rheological mode and seismic cyclic activation that appeared to have been simultaneous with more viscous, probably aseismic, deformational modes (ductile creep). This transitional rheological context, which was contemporary to gold deposition, is supported by field evidence including brittle–ductile auriferous faults (with low- mylonitization) as well as cross-cutting relationships between pure tensile gold-bearing veins and ductile features (folding on various scales and boudinage). The cyclic deformational mode is convincingly indicated by the heterogeneous distribution of fluid inclusions in quartz veins [7] and the ubiquitous occurrence of crack-and-seal textured veins. The successive mineral ribbons in these veins reflect stress and fluid pressure fluctuations as well as repeated opening and precipitation/closure stages induced by seismic cycles [8, 48].
2.2. Relationship to Mineralization
The CLLF contains many gold deposits distributed mostly in clusters that are historical mining districts, most notably the Val-d’Or, Malartic, and Kerr Addison districts [5]. Spatially, they are located next to the CLLF, but not directly upon the fault, and correlate with the major structural segments as described by Bedeaux et al. [41]. Smaller gold deposits also occur along the fault but remain marginal compared to the tonnage found in the clusters (Figure 1(b)). Rafini [30] pointed out in his metallogenic synthesis of the CLLF that these clusters regroup gold deposits having similar characteristics in regard to sulfides abundance and type, alteration, and vein paragenesis as well as orebody shape (see [41] for the English version). This typology of deposits showed the existence of “metallogenic fields” along the CLLF. Several deposits highlight the various parameters, originating either locally or distally, that influenced the deposits. These parameters include fluid origin and chemistry [49, 50], contribution of magmatism [51], host rock geochemistry [52], pervasive primary permeability, and the local structural framework (fracture-related permeability, rheological contrasts, and dilation [49]).
Great effort has been made to place a date on gold mineralization, especially the timing of deformation. At a district scale, gold mineralization appears to have occurred over several events [53–56]. Hydrothermalism was globally active over a long period (tens of millions of years), exceeding the duration of deformation [57]. Still, several authors suggest that a dominant gold deposition event is associated with the late NW–SE shortening event [45, 49, 58–61] expressed predominantly as a dextral strike-slip.
3. Theoretical Background and Methodology
Several models were created to assess the impact of the subsidiary fault framework surrounding the CLLF, the preexisting structural fabrics such as the ubiquitous regional penetrative foliation, the rheological contrasts between lithological units, and the doublet aftershocks.
3.1. Segmentation
During an earthquake, the propagation of slip from the epicentre throughout the rupturing fault window is extremely sensitive to the geometry of the fault [17]. Slip is confined to major segments having consistent geometry and it diminishes drastically approaching the tips [62]. Major segments can also be identified using the slip distributions of modern active faults [17]. Tips of segments occur as either soft links (i.e., gaps or step-overs) or hard links (changes in orientation of fault trace, connecting faults, ramps). During the seismic activation of a fault, these terminations can act as barriers—or rupture arrest points [17, 63]—especially if the termination is a misoriented feature or denotes an abrupt change in geometry. A strong rupture arrest point will therefore recurrently restrain slip propagation and influence static stress distribution and the aftershock location. This will also focus hydrothermal pathways over specific damage zones allowing iterative gold endowment and the formation of deposit clusters. At the scale of the fault, the numerous gold deposit clusters will be related to the repeated failure of several segments [26].
3.2. Static Stress Changes and the Coulomb Criterion
Following a major rupture, the distribution of elastic strain is drastically modified around the displaced fault segment. Corollary changes in stress fields bring some regions closer to a critical failure state, which become the loci of coseismic activity, episodic hydrothermal pathways, and discharge zones [15, 26], while other regions are seismically silent. Geomechanical numeric codes permit simulating, at any point in the model, these static stress changes, induced by sudden displacement, onto a given segment of a main fault and also permit calculating the associated changes in the Coulomb failure stress (e.g., [20, 21]).
The Coulomb-Navier failure criterion is a common empirical measure of the static frictional strength of a fracture. The Coulomb failure stress (CFS) is given bywhere is the change in CFS, is the change in shear stress, is the change in normal stress, is the fluid pressure, and is the friction coefficient. We use the convention that compression stresses are positive. It is a common approximation that, after a seismic rupture, the fluid pressure at any nearby location is not connected to far-field conditions, but rather it is a fraction of normal stress [23, 64]. This theory is referred to as the undrained crust hypothesis, which essentially relies on the fact that static stress transfer is much faster than pore pressure diffusion. This is expressed through the empirical Skempton coefficient :where (2) can now be written:where The Skempton coefficient tends toward 1 when the angle between sigma 1 and the fracture network decreases [65]. The value of μ′ is expected to range from 0.2 to 0.8 [20, 21, 26, 66, 67].
3.3. Doublet
Major tremors are followed by numerous aftershocks that decrease in intensity over time. Aftershocks having magnitude roughly similar to the mainshock are designated as doublets. The segment triggered by this aftershock will in turn influence the local stress field. Thus, the Coulomb stress field can be considerably modified in areas situated far from the mainshock location. An earthquake can be considered a doublet aftershock if [34, 68, 69] the segment experiencing the aftershock is located inside the damage zone related to the mainshock, the two earthquakes have similar magnitude, they have close spatial relationships, and both have close temporal relationships.
3.4. The Numerical Code
The geomechanical simulations were carried out using the UDEC (Universal Distinct Element Code) [70], a finite element code that has the singular ability to interactively simulate internal block deformations as well as displacement along discontinuities that separate these blocks. Although it was developed for engineering purposes, its application to geological problems has proved successful [71, 72]. The 2D model is divided into blocks separated by rheological boundaries that allow the representation of regions having various mechanical properties and discontinuities such as faults and contacts. Internal boundary initial conditions are displacement imposed on the failing fault segment, while far-field boundaries are free. In other words, these conditions are strain-only and there are no stresses applied to the system. The influence of the regional stress state hence is not considered, for the purpose of investigating the coseismic-only stresses and elastic strain, that is, exclusively those propagating from the rupturing segment. Implementing any regional stress field indeed is not relevant to coseismic stress propagation modelling at this stage: either the regional stress field would be presumed initially homogeneous in the entire domain, which implies no relative impact on the results, or the time-dependent and heterogeneous modifications of the stress tensor during the seismic loading-discharge cycle in the vicinity of the fault are modelled, which is a nonvaluable increase of complexity and uncertainty to the scope of this study.
Static stress modelling is a relatively straightforward method where stress and strain are calculated at any point from the derivative of the displacement field converted into the elastic strain field ,and then into the stress field by means of the elastic constitutive law (Hooke’s laws): where and are the Lamé coefficient (Pa) and the shear modulus (Pa), respectively. The symbol represents the Kronecker delta, which is 0 or 1 if . Finally, the resulting stress field is converted into forces applied at each node of the spatial discretization, from which the next-step displacement field is calculated using Newton’s laws.
Finally, the assumption of elastic-only rheological mode into the blocks is a common and consensual first-order approximation [15, 26], highly suitable for the purposes of this modelling. Implementing elastoplastic deformational mode would increase complexity and require additional, weakly constrained, rheological parameters, inducing a barely calibrated and incidentally low (from our testing) influence on the results and then increasing uncertainty in a nonvaluable fashion.
3.5. Geomechanical Parameters
Coulomb stress changes during the failure of the CLLF were modelled in several stages, from simple to more complex and multiparameter model configurations by incrementally including the various geomechanical parameters (Table 1). The simplest model includes only the CLLF (model 1). A structural framework is then added (model 2.a), followed by the integration of failure planes having variable orientations (model 2.b). The influence of potential doublets is finally considered (models 3 and 4). At each stage, the vertical, east-striking Joannes Segment is triggered during the mainshock. This segment was chosen because of its homogeneous geometry and its optimal orientation with respect to the late NW–SE shortening deformation event. Moreover, the segment is clearly delimited by the Davidson Fault to the west and the abrupt bend of the CLLF to the east. Both features constitute well defined and highly plausible rupture arrest points. The Joannes Segment is also located between the Rouyn and Malartic districts, two well-documented areas after decades of exploration. This 41 km long segment is expected to experience earthquakes having 6.6 to 7.0 magnitude associated with 1.4 m maximum lateral offset based on the log-normal empirical relationships established between the earthquake magnitude, the length of the ruptured segment, and the offset (Figure 3, [31, 32, 73, 74]).
In each model, blocks are assigned geomechanical properties corresponding to their dominant lithology (sedimentary, volcanic, or granitoid rocks). These are based on the abundant experimental data available in the literature [75–78]. Block boundary properties are also determined according to their structural nature: major faults, second-order faults, third-order faults, or simple contacts (Figure 2). Such a hierarchical assignment of properties is designed to depict the actual relative importance of every contact in the deformational framework, based on an extensive knowledge of the field area obtained over several decades of observations in southern Abitibi [40–45, 79, 80].
The CFS may be calculated for any point from the projection of the local stress tensor along the failure plane of interest. Consequently, a critical parameter is the orientation of failure planes, which relies on various structural assumptions. A common assumption is to consider primary failing planes as being oriented optimally at every point of the model to maximize Coulomb stress [20, 22, 23]. This means setting the angle of as constant and equal to 45–0.5 atan (μ) (i.e., to 29.5° for ). This method can be suitable if the model includes isotropic rocks having no preferential failure plane and if regional stress tensors are well defined. In the case of the CLLF, the late strike-slip event is related to a general NW–SE shortening event where the exact orientation is not defined. This means that the optimal failure plane orientation can only be estimated, and this estimate can range from N090° to N120° with respect to the direction that can range from N300° to N330°. An alternative assumption is that preexisting, low-strength fabrics preferentially undergo failures. Around the CLLF, such fabric corresponds to the ubiquitous and highly penetrative regional foliation, which is mostly east-striking and globally parallel to bedding and faults [36, 41]. Its similarity to the optimal orientation makes it a reasonable failure plane alternative, and this is supported by field evidence showing reactivation during late strike-slip events [41, 49]. Moreover, the existing surface geology in southern Abitibi corresponds to a crustal level close to the one experiencing high-magnitude earthquakes when the CLLF was active. In an initial stage, failure planes were set as east-striking by default. In a second stage, the refinements to the model obtained through implementing foliation fabric orientations at any point were investigated through the use of an interpolation of the numerous available foliation measurements (Figure 2). For calculating the CFS, a value of 0.4 was used for μ′ based on previous modelling attempts [20, 26, 29].
The final investigated parameter was the influence of doublet aftershocks. In the vicinity of the CLLF, the Rivière-Héva and Lapause second-order faults are high-potential doublet segments using empirical criteria based on present-day seismic observations: location in a zone of positive CFS change and close proximity to and similar segment length as the Joannes Segment experiencing the mainshock. The first segment considered to be able to produce a doublet aftershock is associated with the Lapause Fault, which marks the boundary between sedimentary and volcanic rocks in the eastern part of the CLLF (Figure 2). The 24 km long segment trends at 120° and is delimited by an abrupt change in azimuth in the northwest; it intersects the third-order Parfouru Fault in the southeast. Coseismic failure along this segment is expected to yield 1.3 m of dextral displacement based on its length (Table 1). The other segment is related to the Rivière-Héva second-order fault, a recently discovered structure in the volcanic rocks of the Malartic Group, which has significant importance in the local stratigraphy [79, 80]. The selected segment strikes 115° to southeast and is bounded on each tip by an abrupt change in azimuth to the east. Its length of 29 km corresponds to dextral displacement of 1.4 m for one earthquake.
4. Results
Static stress change mapping around the CLLF was done by the testing of seven models (Table 1 and Figure 2). Three models involve failure on the Joannes Segment of the CLLF with extra parameters added. The first model uses only the CLLF. The second model adds the geological and structural frameworks of southern Abitibi, that is, additional faults, contacts, and major lithological assemblages. The third builds on the second model with the addition of variable orientations of the failure planes for CFS calculations. Potential doublets along the Lapause and Rivière-Héva Faults were initially tested using the geological and structural frameworks of Abitibi and then again tested with the addition of variable failure plane orientations. Models involving doublets are combined with the CFS calculations from the Joannes Segment failure.
4.1. Cadillac-Larder Lake Fault Failure (Model 1)
In the generic model 1.a, including only the CLLF, lobate areas representing positive and negative change in Coulomb stress are distributed in a symmetric pattern, similar to the results obtained by King et al. [20]. Positive lobes, corresponding to predicted damage zones, are distributed at the tips of the fault failure plane. Negative lobes are located on each side of the Joannes Segment (Figure 4(a)). The eastern positive lobe is well developed with an ESE trend, although the W-trending western lobe is less developed.
(a)
(b)
(c)
Model 1.b (Figure 4(b)) incorporates southern Abitibi geological and structural frameworks. Lobes are generally smaller than those within the generic model. The main difference with the first model is that lobe boundaries are heavily influenced by contacts and faults. In addition, lateral positive lobes are absent. The western positive lobe is delimited by the Davidson Fault to the east and covers a larger area than in the generic model. It also extends farther west along the CLLF.
The final model (model 1.c) to represent the main rupture along the Joannes Segment displays the most important variations from the classic geometry of Coulomb stress change only when the preexisting low-strength fabrics are included as the primary potentially failing planes (Figure 4(c)). Taking into consideration the variable orientation of the failure planes induces a locally important increase or decrease of Coulomb stress change. The eastern lobe strikes SE and shows a restricted width and length. It lies close to the CLLF next to the tip of the Joannes Segment. The western damage zone is larger than those within the previous models and is well developed next to the intersection between the Davidson Fault and the CLLF. In contrast with other models, it extends to the south of the CLLF instead of being mostly confined to the north of the fault.
4.2. Doublet-Lapause Fault (Model 2)
Model 2.a includes the geological and structural framework and generates numerous lateral positive lobes in addition to lobes located at the tip of the ruptured segment (Figure 5(a)). High positive CFS occurs between the Joannes and Lapause Segments. The eastern lobe covers a larger area and is east-striking. Lateral lobes are north-striking and are also associated with southeast-striking faults proximal to the Lapause Fault.
(a)
(b)
Model 2.b involves variable failure plane orientations (not shown). Positive CFS changes are significantly affected by variations in the orientation of the failure planes and are mostly confined to the tip of the segment. The eastern lobe strikes WNW and is developed along the east-striking faults to the north of the CLLF. The western lobe is nearly absent and is confined to the SE-trending section of the CLLF.
4.3. Doublet-Rivière-Héva Fault (Model 3)
In model 3.a (Figure 5(b)), only the lobes at the tip of the Rivière-Héva segment are developed. The western lobe is east-striking and delimited to the south of the Rivière-Héva Fault and reaches the tip of the Joannes Segment. The eastern lobe closely follows the east-striking splay fault of the Rivière-Héva Fault.
In the final model, the positive lobe between the Joannes and Rivière-Héva Segments is poorly developed and is confined to the SE-trending portion of the CLLF (not shown). The eastern lobe is much more developed and covers the area between both splay faults of the Rivière-Héva Fault. A N-trending lateral lobe is also well developed from the western tip.
5. Discussion
In this section, we review the validity of the method and compare results of models among each other. We also discuss implications for gold metallogeny and exploration in southern Abitibi.
5.1. Validation of Methods
5.1.1. Frequency-Magnitude Diagrams
Following ruptures in the main fault segment, lower-scale displacement and lower-magnitude seismicity are expected to occur along major adjoining and surrounding fault segments [21]. Simulations integrating the structural framework partially allow this slip in accordance with the faults’ relative importance as determined from field observations. This suggests implicitly that failure conditions are reached on these faults. Displacement fields calculated on contacts in the entire domain following the main rupture are converted to frequency-magnitude profiles by means of classical empirical relationships [31, 32]. The profiles obtained from simulations are then compared to actual profiles to confirm the mechanical validity of the models (Figure 6). Frequency-magnitude profiles are commonly represented by the Gutenberg–Richter relationships:where is the number of rupture events with moment magnitude greater than and and are constants [29, 33, 81]. This power law behaviour was shown to be related to elastic interactions within a heterogeneous medium (e.g., [82]), with increasing with heterogeneity. The value of has been the focus of extensive interest, primarily for seismic hazard prediction, as it was recognized early as a highly sensitive parameter in crustal rheology. At the scale of the fault system, it is expected to vary between 0.6 and 1.2 depending on various factors including temporality relative to the seismic cycle (e.g., [29, 83–85]). Other controlling factors are more context-related and are of greater interest in validating our models: the stress regime such as type of fault, depth [86], and the slipping mode—notably the creeping mode which increases [33]. We reviewed values recorded on large transcrustal active strike-slip faults in oblique convergent margins, which may be considered as modern analogues of the Archean CLLF deformations modelled in our study: the San Andreas Fault ( and [87, 88], resp.), the North Anatolian Fault ( [88]), the Central Japanese Median Tectonic Line ( [87]), and the Sagaing Fault area, Myanmar ( [85]). Perceptible variations are observed partly due to methodological bias, but recorded values consistently trend below unity, mostly ranging from 0.5 to 0.9. Simulated profiles display values ranging from approximately 0.79 to 0.84 (Figure 6), which are similar to actual trends, indicating that our models reliably reproduce natural frequency-magnitude distributions for the targeted context. Although this mechanical validation attempt is conclusive, it overrides the temporality of values when aftershock sequences are compared to historical records. This is an issue requiring a verification in future work.
5.1.2. Correlation with Known Gold Deposits
Gold deposits in southern Abitibi are clustered near the CLLF, forming several goldfields as illustrated in Figure 1. Damage zones are identified as a critical control in the generation of gold deposits as they allow sustained fluid flow while allowing precipitation of gold-related minerals [15, 24, 29]. Previous static stress mapping studies in the Archean goldfields in Australia have shown that the clustered distribution of gold deposits along major fault systems can be correlated to positive CFS changes, a proxy for damage zones, following failure along specific fault segments (e.g., [9, 26]). In southern Abitibi, exploration for gold deposits has been carried out to various degrees using geological features such as major faults as criteria for gold prospecting. This yielded abundant documentation, sampling, and drill hole surveys as shown in Figure 7. This figure displays fewer exploration studies in the Pontiac Subprovince due to a thick glacial overburden and exploration strategies traditionally favouring volcanic rocks that are largely absent south of the CLLF. Only a few, less than a kilometre wide areas have been explored in detail; all correspond to gold deposits and gold mines (i.e., the Joanna gold deposit, Sigma-Lamaque, Canadian Malartic, Lapa, and Laronde gold mines). Overall, the known distribution of gold north of the CLLF does not correlate with the intensity of exploration studies. Moreover, these are well-studied regions, as there are no areas lacking in data, as shown in Figure 7. As such, the present-day spatial distribution of gold deposits and tonnage can be realistically considered as unbiased indicators of fossil auriferous hydrothermalism north of the CLLF and in its nearby southern vicinity. Thus, these are reliable validation variables for the models.
(a)
(b)
Models can be compared in regard to their ability to match the distribution of existing gold deposits and tonnage. For a model to be deemed a good fit, it is expected that damage zones correlate with known gold occurrence. The detailed results of this validation analysis are given in Table 2. Positive CFS changes are considered significant when above 30 kPa [20, 21].
The Rouyn goldfield is only reproduced by damage zones induced by rupturing of the Joannes Segment. In all simulations involving failure of the Joannes Segment, 73% of the Rouyn goldfield deposits are successfully predicted as they are correlated to significant positive CFS change. However, models 1.b and 1.c offer a better tonnage correlation, accounting for 48% of the existing tonnage. The Bousquet goldfield is poorly explained by the Joannes Segment failure, with only two deposits (11%) located west of the goldfield covered by the damage zones. Such results were expected, as the triggered segment (containing the Bousquet district) experiences an important stress drop following failure [20]. In addition, the Bousquet goldfield is located quite distant from doublet segments.
The Malartic goldfield shows a greater variability between models. The model that includes only the CLLF produces a correlation with 90% of existing deposits. On the other hand, the model incorporating the Abitibi structural framework and the doublet aftershock on the Lapause Fault can account for only 67% of existing gold deposits, suggesting these factors are not significantly involved in gold genesis in this area. In all cases, tonnage correlations to damage zones in the Malartic goldfield are markedly low. This is attributed to a bias due to the important proportion of gold contributed by the Canadian Malartic mine, a dominant intrusion-related deposit accounting for 67.5% of the total Malartic goldfield tonnage. Smaller deposits located over a wide area along the CLLF still occur in positive lobes.
In all three models involving only the Joannes Segment, the Bourlamaque and Val-d’Or goldfields remain largely unexplained by static stress changes. Moreover, static stress is expected to be modified only 10–20 km away from the failing segment [20, 26]. This means that gold deposits in the Val-d’Or and Bourlamaque districts cannot be explained solely by failure of the Joannes Segment as they are too far away from this segment. In the case of the Val-d’Or segment, the doublet failure of the Lapause segment matches 38% of existing gold deposits and 95% of the goldfield tonnage in both models 3.a and 3.b, revealing the importance of this segment for explaining the eastern deposits along the CLLF. More precisely, this doublet rupture corroborates Val-d’Or goldfield’s most important deposits such as Sigma-Lamaque, Goldex, and Kiena. It is worth noting that the cluster of deposits located near the recently developed Marban deposit is also remarkably well predicted. Doublet rupture along the Rivière-Héva Fault yields damage zones containing 40% of gold deposits in the Bourlamaque district and 50% when the variable orientation of regional foliation is integrated (model 3.b). These results illustrate how the doublet seismic activation of the Rivière-Héva Fault—a subsidiary of the CLLF—extends the influence of the CLLF laterally and permits gold occurrences in the Bourlamaque district despite being located quite distant from the CLLF.
In summary, the scenario involving the sole failure of the Joannes Segment produces a good correlation with deposits in the Malartic and Rouyn districts. However, a combined model that simulates successive rupturing of the Cadillac-Larder Lake and Lapause and Rivière-Héva Faults produces the best results for explaining the occurrence of gold in regard to the number of deposits and related tonnage. During a deformational event, faults will be active over millions of years, involving multiple repetitions of displacement during rupture events. Each of these events will trigger neighbouring major and second-order fault segments, including doublet aftershock ruptures on the Lapause and Rivière-Héva Faults. Over time, fluid flow will be focused along the same pathways, progressively developing and endowing gold deposits. In southern Abitibi, repeated doublet triggering of second-order faults by mainshocks on the Joannes Segment allows gold deposition to occur in areas farther away from the CLLF. It is thought to be possible that some gold deposits were developed by this combination of several failing segments repeated over time. These models demonstrate the significant role of the regional geological framework and orientation of the regional foliation on the spatial distribution of the positive lobes and how doublet aftershocks play a significant role in explaining the location and endowment of existing goldfields.
5.2. Hydrothermal Fields Validation
In a continued effort to constrain the models’ validity, this section explores the spatial correlation between hydrothermal fields derived from geochemical and mineralogical data and the modelled damage zones expected to be triggered due to coseismic permeability enhancement. A critical control for gold precipitation is the destabilization triggered by wall-rock chemical interactions [89, 90] and by fluid pressure changes induced by fault cycling [6, 29], as explained in Section 5.5. Mesothermal gold-bearing fluids are in a geochemical disequilibrium with the host rocks [10, 24, 28], resulting in fluid-rock interactions that primarily include transfers of H, C, and S elements from the fluid to the wall rock by replacement or mineral precipitation. These exchanges cause substantial hydration, carbonation, and, to a lesser extent, pyritization. These common reactions, along with silicification, albitization, and potassic alterations, mobilize major cations (Ca2+, Mg2+, Fe2+, Si4+, Na+, and K+), which are derived either from the in situ destabilization of wall-rock minerals (e.g., [89, 91]) or from external metasomatic input supplied by the fluids. Finite-state absolute gains and losses therefore provide, to some degree, records of the composition of the hydrothermal fluids that percolated over time. Thus, attempts to identify different signatures of fluids and delineate distinct hydrothermal fields are made possible by calculating multielement mass balances on a regional scale from existing lithogeochemical databases.
Absolute gains or losses of Na and K have been found by Bigot [92] to bear the highest potential for ore classification in southern Abitibi, since most mesothermal gold deposits exhibit significant enrichments in either Na or K, although rarely both. Carbonatization and silicification, for instance, are ubiquitous in the CLLF’s surroundings. Na-K metasomatism is seemingly more limited in space and time and is related to auriferous hydrothermalism. Rafini [30] mapped Na and K mass balances in the Malartic and Val-d’Or districts using a compilation of 3488 whole-rock analyses of mafic to felsic igneous rock samples from government (SIGEOM), academic, and industry databases. Mass balance calculations use the precursor-modelling method [93], in which major cation compositions of fresh rocks are predicted for each sample using the ratios of immobile elements Al, Ti, Y, and Zr based on relationships established by neural networks on a large worldwide reference database of unaltered rocks. This method is analogous to isocon in that immobile elements permit predicting cations gains or losses, yet the precursor (unaltered) rock composition is modelled rather than measured (sampling unaltered rocks is not required). The results (Figure 8) suggest logically organized zones that evoke distinct hydrothermal fields, yet their delineation is not straightforward. This is partly due to a highly heterogeneous sampling coverage; a given cell value may correspond either to a single sample or to a predominant trend derived from fifty samples. Moreover, strong local-scale variations are visible, which plausibly reflect the superposed footprints of several alteration events associated with the long-lived CLLF hydrothermal system. Notwithstanding these reservations, such analysis of the finite-state metasomatic regional footprint provides a valuable overview of the predominant hydrothermal events. Eastern regions, where structural fabrics (stratification and schistosity) are east-striking, record large-scale Na-K losses along with very strong local-scale K-gain zones; these patterns have been interpreted as representing a volcanogenic footprint [30]. It is interesting that this metasomatic pattern changes drastically west of the “axial zone” where structural fabrics shift from E- to ESE-striking. Our numerical simulations predict that the influence of the Joannes Segment ruptures ends approximately in this “axial zone” (Figure 4). This suggests that fluid flow triggered by persistent rupture arresting at the Lapa bend would not reach the regions east of the “axial zone,” hence predicting a geographical break in the metasomatic record that is consistent with these observations. This also echoes the fact that Val-d’Or and Bourlamaque goldfields, which are east of the “axial zone,” show apparently distinctive structural characteristics, with the predominance of subhorizontal tensile veins associated with E–W auriferous faults (discussed in Section 5.3). West of the “axial zone,” the presence of distinct ESE-oriented fields north and south of the Lapause Fault is noticeable: the latter exhibits a drastic K-gain corroborated by recurrent microcline occurrences in mineralized rock samples as well as a weaker and less ubiquitous Na-gain. In contrast, the former shows a large K-loss associated with a strong and widespread Na-gain that superposes markedly on deposits where a significant albitization is reported in mineralized veins. These hydrothermal fields match the geometry of the damage zones obtained from our doublet rupture simulations along the Lapause Fault (models 2.a and 2.b, Figure 5(a)). These models predict a dissociation in space and time between hydrothermalisms occurring north (deposits proximal to the Marban mine) and south (deposits proximal to Malartic mine) of the Lapause Fault. These models reflect the occurrence of distinct metasomatic signatures along these two areas.
(a)
(b)
5.3. Comparison with Other Ore-Targeting Models and Limits
This study focused on modelling CLLF seismic activity using existing geological knowledge. In doing so, several assumptions were made, which limited our results although they offer potential future research topics. Modelling late strike-slip events ensures that fault geometries remained largely unmodified at a regional scale. However, the model remains largely dependent on our current knowledge of regional geology, most notably the traces of faults and their relative importance. In particular, the geometry of southern Abitibi at depth remains largely speculative and will remain an obstacle for precise and representative 3D modelling of fluid circulation and its influence on potential discharge areas. Another critical limitation is the orientation of failure planes. Even if foliation planes surrounding the CLLF are mostly east-striking, variations can still occur at a local scale along shear zones, contacts, and disrupted foliations, which are not accounted for in our interpolation.
Several deposits in the Val-d’Or and Bourlamaque districts (including Sigma-Lamaque), correlating with predicted damage zones, have previously been interpreted as formed during reverse movement on E–W faults associated with N–S shortening [94–96]. However, faults recorded in these deposits mostly display oblique reverse movements that are compatible with NW–SE shortening [97, 98], meaning that part of the mineralization in these districts can be related to the late dextral strike-slip deformational event. Moreover, damage zones in Val-d’Or and Bourlamaque districts are triggered by the Rivière-Héva and Lapause Faults, both of which are southeast-striking and steeply dipping. It is thus possible that these two faults were also triggered as dextral faults during N–S shortening.
The model that integrates the regional framework, failure along the main fault, and seismic doublet on second-order faults explains very well the distribution of gold deposits in the Rouyn, Malartic, Bourlamaque, and Val-d’Or districts. However, the Bousquet district remains largely unrelated to positive CFS changes in a dextral strike-slip fault system. Gold deposits in this district are distributed more homogeneously along the CLLF, which obviously differ from the other districts in which most of the deposits are in the surroundings of the CLLF. This likely reflects another form of structural control in the Bousquet district, which comforts the hypothesis that the Joannes Segment is the rupturing master segment exerting a prevailing control on the formation of goldfields at a regional scale. This also corroborates our results in that no damage zone is generated along the rupturing segment itself, which explains why no deposit formed around the fault in the Joannes Segment.
5.4. Coseismic Release of Deep-Seated CO2 Reservoirs
In seismically triggered hydrothermal systems, the episodic enhancement of permeability along crustal-scale fault drains creates a transient hydraulic continuity between overpressured, sublithostatic deep reservoirs and possibly hydrostatic, supracrustal fluids (e.g., [26, 99]) (Figure 9). Supracrustal levels are the locations of orebody formation. In the CLLF area, deep-rooted connections between faults have been tentatively documented [52, 98]. Failure can thus reach deep pressurized reservoirs and release fluid flows through the entire fault system.
Carbon-dioxide-rich metamorphic fluids have deep sources originating from crustal devolatilization (likely from tectonically buried carbonate sediments) during prograde metamorphism [100, 101]. The high CO2 content in gold-bearing fluids, commonly from 5 to 20 mol%, has been reported extensively from fluid inclusion analyses [2, 101–103], including into deposits related to the CLLF [48, 104, 105]. These fluids are associated either with late orogenic midcrust reheating or with subduction-related processes and subsequent fluid entrapment in the mantle wedge [5]. In the region of the CLLF, it is likely that the abundant CO2 inflow, reflected by ubiquitous carbonate alterations (e.g., [40, 98, 106, 107]), is supplied by vertical fluid transfers from subamphibolite devolatilized Pontiac sediments enrooted at the CLLF’s footwall (Figure 9), as suggested by the LITHOPROBE profiles [108]. It is worth noting, however, that both midcrust and mantle fluid source models are remarkably well supported by present-day observations of crustal fluid transfers following major earthquakes in active convergent zones. Husen and Kissling [109] and Nippress and Rietbrock [110] reported the upwelling of a fluid mass and a fluid pressure front, respectively, through the entire crust thickness following major thrust-fault earthquakes in Chile. Miller et al. [111], Antonioli et al. [112], and Terakawa et al. [113] observed that the seismic activation of major faults in Switzerland and Italy consistently triggered the release of entrapped CO2-rich overpressured fluids from the base of the seismogenic crust to the upper levels; the poroelastic diffusion of these fronts is considered to be a first-order control on aftershock spatiotemporal sequences. Such cyclic CO2 pulses agree remarkably well with the genetic model of mesothermal gold deposits, in which veins develop from pulsatile infilling of CO2-rich multiphase fluids, as well as with the extensive carbonatization associated with fossil hydrothermal gold systems.
5.5. Multiphase Fluid Flow and Metal Deposition
Phase separation likely takes place during hydrothermal coseismic triggering and the transient upward migration of metamorphic fluids. Phase separation has been proposed as a first-order mechanism in mineralization processes related to hydrothermalism [114, 115], primarily magmatic-hydrothermal systems [116], and mesothermal gold [24, 117, 118]. During the seismic cycle, deep-seated fluids undergo drastic adiabatic depressurization associated with its migration through brittle crustal levels, which triggers phase separation and the partitioning of the CO2-rich vapor phase. Experimental works demonstrated that the fluid’s CO2 content exerts a significant control on the temperature and pressure of phase separation [115, 119], which may hence play a governing role in the depth of gold precipitation. These authors indeed highlighted that CO2-free fluids would undergo phase separation at higher crustal levels than that commonly reported for mesothermal gold deposits. Pokrovski et al. [119] concluded from thermodynamic constraints that Au is partitioned into the vapor phase in the presence of reduced sulfur in acidic-to-neutral conditions, a possible enrichment process associated with phase separation [115]. Zezin et al. [120] demonstrated that the gold solubility in the vapor phase is due to hydrated sulfide complexes, whose stability is closely dependent on H2O partial pressure. In contrast, Hodkiewicz et al. [117] submitted a possible gold precipitation trigger related to phase separation where the partitioning of reduced gas in the vapor phase modifies the oxidation state of the liquid, which destabilized bisulfide complexes and hence gold solubility. Additionally, as a possible gold precipitation trigger, the partitioning of reduced sulfur in the vapor during phase separation would induce a drop in sulfur fugacity which would also destabilize Au bisulfide species from the liquid. Other recent works have documented the decrease in noble metals solubility as a direct consequence of fluid pressure drop and phase separation [24, 121].
In this system, the vaporization of fluids induced by pressure drop in damage zones is a plausible process exerting control on Au solubility and gold deposition, possibly involving changes in oxygen fugacity as explained above. It is worth noting that other parameters have the ability to trigger Au deposition: any change in pH, oxygen, sulfur, or hydrogen fugacity related to wall-rock chemical interactions [89, 90] would tend to destabilize gold complexes either from the liquid or from the vapor phase (e.g., [5, 101, 117]). Although the detailed distinct implication of these various processes and their interrelationships in the formation of mesothermal gold ore still give rise to discussions, it is likely that several processes rather than a sole one are involved at a district scale [30] and even at a deposit scale.
5.6. Poroelastic Modelling
This study focused on the role of coseismic static stress changes in the genesis of hydrothermal gold deposition by hypothesizing a role for undrained crust conditions. However, another aspect is the influence of the transient diffusion of fluid pressure triggered by earthquakes. Indeed, seismically triggered fluid circulation is expected to result in marked pore pressure variations that also modify the stress field and the geometry of potential fluid discharge zones. This presumably also triggers aftershocks in areas not experiencing sufficient positive CFS change [23, 111, 122]. Such conceptual statements about the coseismic interactions between an overpressured diffusive fluid and the triggering of aftershocks are remarkably well corroborated by recent work in seismology. The observation of present-day seismic events from various geodynamic environments demonstrates that the outward propagation of the triggering front from the main rupture zone is governed by the transient diffusion of overpressured fluid [111, 112, 123–126]. This was strongly suggested by records of the spatiotemporal propagation of triggering fronts in several present-day seismic events (Umbria Earthquake, Italy, 1999; Vogtland Swarm, Eastern Europe, 2000; Northern Apennines Sequence, 2000; Antofagasta Earthquake, Chile, 1995). The distance travelled by these fronts exhibits a square root relationship to time, reflecting a classical diffusive law [123]. This indicates that the spatial distribution of aftershocks is at least partially controlled by midcrustal fluid pressure diffusion. A consequence is that the undrained crust hypothesis, related to the use of Skempton’s parameter (see (2)), may not hold on the time scale of aftershock occurrences (weeks, months, etc.). Static stress transfers are practically instantaneous whereas the fluid pressure relaxation extends over months and clearly reproduces the observed time decay of aftershocks triggering. As a result, the maintaining of enhanced permeability—and related hydrothermal fluid transfers—near the main rupture zone may actually be the result of coupled effects of purely static stress transfers and stresses related to the transient diffusion of fluid pressure. This is the proper definition of the poroelastic approach [23]. Based on the Landers 1992 earthquake dataset, these authors demonstrated that poroelastic models provided stress mapping having a slightly better agreement with real aftershock locations than did the static stress models (87% versus 65% of epicentres were located in positive Coulomb stress change zones). However, due to the configuration of the poroelastic model and the close interactions between both stress propagation mechanisms, the two models displayed strong similarities. In conclusion, even though static stress models may explain, to a large extent, the spatial distribution of aftershocks [127], further study should consider incorporating poroelastic relaxation into the modelling of coseismic stress mapping in an attempt to better reproduce fossil hydrothermal fields.
6. Conclusion
This study proposes a quantitative framework where repeated seismic ruptures of a crustal window recurrently trigger fluid ascent through temporarily permeable damage zones via the breaching of the fluid pressure differential between the deeper and the upper crustal levels (Figure 9). These zones are located all around the ruptured surface (that extends over the entire seismogenic crust width for earthquakes , [32]). Iterative ore deposition is produced from repeated fluid discharges into the affected areas and is also related to persistent rupturing of the same fault window. This occurs when the fault shows a strong structural segmentation, leading to rupture arrest points on both extremities that remain constant over multiple seismic events. In the Joannes Segment, such robust rupture arrest points are the drastic inflexion of the CLLF around the Lapa deposit (east tip) and the early offset on the Davidson Fault (west tip). Our models markedly suggest that these structural barriers are responsible for the formation of goldfields in the Malartic and Rouyn districts.
The influence of various mechanical and rheological factors on the spatial distribution of coseismic damage zones was also investigated. The role of these factors in coseismic gold mineralizing processes was evaluated through the ability of the models to corroborate known deposit locations. The models tested in this paper explained the distribution of many gold deposits spatially associated with the CLLF: 60% of all gold deposits, including 90% and 73% of the Malartic and Rouyn goldfield deposits, respectively, and 48% of the total gold tonnage, including 95% and 71% of Val-d’Or and Bourlamaque goldfield tonnage, respectively. Implementing geological and structural frameworks improves the models’ correlation with known gold deposits producing a more detailed delineation of zones of hydrothermalism triggering. Doublet seismic ruptures along second-order faults are shown to be important factors in the genesis of goldfields at distances away from the CLLF. It is shown that the combined seismic events along the CLLF and, alternately, the Lapause and Rivière-Héva Faults produce extended damage zones, which explain gold deposition in the Bourlamaque, North-Malartic, and Val-d’Or areas, including the Sigma-Lamaque, Goldex, and Marban major deposits.
Identified hydrothermalism triggering zones, as predicted by the various combinations of main or doublet rupturing segments, show compelling geographical similarities to hydrothermal fields identified through the analysis of the regional metasomatic footprint using mass balance calculations and mineralogical observations. This provides a first-order, valuable validation of the submitted concepts.
The Malartic and Val-d’Or/Bourlamaque goldfields, that is, located, respectively, west and east of an axial “kink” zone where the structural fabric trends change from ESE–WNW to E–W, relate to distinct hydrothermal systems in time and space. This is suggested by the coseismic damage influence of the Joannes Segment rupture ending up at this axial zone. This gives perspective to the fact that Val-d’Or and Bourlamaque goldfields differ from the Malartic goldfield in terms of bulk metasomatic footprint and structural features.
The models are geomechanically validated by calibrating frequency-magnitude seismic distributions of the simulated aftershock sequences with that of actual datasets from analogous active major strike-slip faults such as San Andreas and the North Anatolian systems.
Finally, a likely source of the seismically migrating fluids is the devolatilization of the Pontiac carbonaceous platform that is enrooted at the CLLF’s footwall in supra-amphibolite crustal levels; this matches recent conclusions by Gaboury [105] and Leventis [128].
The promoted coseismic hydrothermal plumbing system was quantified and validated based only on its superior expression. Further work will benefit from the ongoing deep crustal geological acquisition programs (Metal Earth Project) to investigate this model over extended crustal windows owing to reliable constraints on the deep geology of the Superior craton.
Additional Points
Highlights. (i) The role of large Archean earthquakes in the genesis of mesothermal goldfields is investigated in southern Abitibi. (ii) Coseismic fluid-focusing damage zones induced by ruptures of the Cadillac-Larder Lake Fault (CLLF) are simulated using static stress modelling. (iii) Persistent seismic ruptures along the Joannes Segment exert a large control on the deposition of the Rouyn and Malartic goldfields. (iv) Seismic doublets along subsidiary faults are shown to be prime factors in generating goldfields geographically distant to the CLLF. (v) Simulated damage zones correlate well with the distribution of gold deposits and the regional metasomatic footprint.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This paper forms part of the first author’s Ph.D. thesis carried out at the Université du Québec à Chicoutimi. The work was supported by CONSOREM, the Mineral Research Exploration Consortium, which provided financial and logistical support, and by the Ministère de l’Énergie et des Ressources Naturelles, Québec, Canada.