Chapter 2
2.1 Introduction
Worldwide, numerous attempts have been made at explaining why, in a given region, crustal earthquakes do not occur beneath a certain depth. To be compared with the rheological characteristics of the crust, earthquake focal depths are evaluated in three different fashions in the literature: the depth where a decrease in the distribution of earthquakes occurs (DDec; Meissner and Strehlau, 1982); the depth above which 80% of earthquakes occur (D80%; Smith and Bruhn, 1984; Doser and Kanamori, 1986; Banda and Cloetingh, 1992); and finally the maximum focal depth of earthquakes (named thereafter DMax ; Chen and Molnar, 1983; Kusznir and Park, 1984; Ranalli and Murphy, 1987; Chen, 1988). As described in Chapter 1, CSZ earthquakes occur exclusively in the mid- to upper crust, from the surface to a depth of 30 km, with DDec at about 13 km and D80% at 15 km (Figure 2.1). In the CSZ, DMax is at 29 km; however, only 1% of earthquakes occur between 25 and 29 km depth. We therefore adopt the depth of 25 km (D99%) as a representative value for the lower limit of the seismic crust. This should not be interpreted as a sign that earthquakes between 25 and 30 km are negligible; as a reminder the 25 November 1988 mbLg 6.5 Saguenay earthquake occurred in the Grenville Province about 100 km west of the Charlevoix Seismic Zone, with an unusual focal depth of 29 km (North et al., 1989). The final discussion of this chapter will briefly examine the significance of this earthquake, and of others in the 25 to 30 km range, in light of the Charlevoix modelling.
Rheological modelling can provide constraints to the earthquake depth distribution, from the seismogenic brittle layers at upper to mid crustal depth to the plastic regime in the mid- to lower crust. To model the factors that favour or inhibit the occurrence of earthquakes, lithology, structure, temperature and pressure, pore fluid pressure, stress and strain rate have to be considered (Ranalli, 1991). In this chapter, some important factors in modeling the brittle and ductile parts of the crust are reviewed and applied to the CSZ. Based on these results, the factors controlling the CSZ depth distribution are discussed.
2.2 Temperature
Crustal temperatures can be estimated assuming that conductive heat transfer is predominant in continental intraplate environments such as Eastern Canada. Steady state 1-D heat transfer in a medium with heat production A(z) is used and can be described by a Poisson equation (see e.g. Chapman, 1986)
where T is temperature, z depth, A heat production rate, and k thermal conductivity. By double integration with respect to z, and by assuming qT is the heat flow at the surface, one obtains for a layer of constant heat generation and thermal conductivity
where TT is the temperature at the surface.
If a layer has thickness Deltaz, then, at the bottom of the layer, the temperature TB and the heat flow qB can be expressed in terms of the temperature and heat flow at the top of the layer (TT, qT) and properties (A,k) of the layer
Equations (2.3) and (2.4) can be used to calculate T(z), resetting TT and qT at the top of each new layer with the values TB and qB solved for the bottom of the previous layer. This formula can be used for one-dimensional modelling of the geotherms where heat production is assumed constant within layers. Although one-dimensional modelling represents an oversimplification (Ranalli, 1991), the uneven data distribution for the Grenville Province restricts studies to the one-dimensional case.
There has been considerable debate as to how heat generation is distributed with depth. Another possibility is an exponential decrease of heat generation with depth
where DA0 represents the contribution from radiogenic heat sources in the upper crust, with D being related to the characteristic thickness of the heat generating layer. This approach is based on the known enrichment in K, Th and U of upper crustal rocks. The exponential decrease model is not considered appropriate anymore; heat generation data from the 12 km deep Kola peninsula well support a step-like distribution rather than an exponential distribution (Drury, 1989). Thus, a step-like distribution is assumed in the following sections.
2.2.1 Surface heat flow
The most important factor in geotherm calculations is surface heat flow. At a given site, measurements are subject to experimental errors (generally considered small) and to near-surface perturbations of conductive heat transfer, such as ground water circulation. In continental environments, surface heat flow can be measured in different geological environments such as in deep wells in hard rock, in dry water wells and in lake bottom sediment columns. At a given site, the assigned heat flow value represents an average of values (temperature gradient and conductivity) from different depths. In general, the deeper measurements are the most reliable, being remote from near-surface perturbations caused by water circulation, topography variations or paleoclimatic disturbances. Due to these factors, it may be necessary to assume an uncertainty of at least 10% in heat flow measurements (Drury, 1987). In assigning a heat flow value to a given geological province, many assumptions have to be made: that quality of the data is high; that heat flow is relatively constant over the area (or that the mean and standard deviation of heat flow values provide a fair estimate of the real values); and that all geological settings of the province were properly sampled. In other words, a heat flow value can be considered typical of an area if heat flow measurements are reliable (i.e. the effects of near-surface perturbations are removed), numerous enough to be statistically significant and spatially distributed over most sub-areas.
2.2.2 Heat production
Although the uncertainty in heat generation for a given sample is relatively small, variability of heat generation can be very high for a given rock type (Drury, 1987). Heat production studies provide estimates of various "typical" rock types that can be related to crustal levels. In general, upper crustal rocks have higher heat production values than the lower crust (generally higher than 0.4 µW/m3, and up to more than 5.0 µW/m3). This range is related to the more felsic composition of the upper crust with its higher concentration in radioactive elements. In most studies, upper crustal heat production is determined from values of surface rocks and from geological and geophysical knowledge of the area. In general, lower crustal heat production values are assumed to be in the 0.4 ± 0.3 µW/m3 range, which corresponds to granulite facies rocks (Ashwal et al. 1987).
2.2.3 Thermal conductivity
Although thermal conductivity may vary from sample to sample, it is less variable that heat production and is generally assumed to be constant within a given crustal layer (generally between 2.5 and 3.0 W/m•K). As will be examined in more detail in other sections, small variations in thermal conductivity have little effect on geotherms.
2.2.4 Thickness of layers
Geophysical information, such as seismic and gravity data, can be used to assess the thickness of the crustal layers. For geotherm calculations, typical values of heat generation and thermal conductivity can be associated with the assumed geological units. Major uncertainties arise when no geophysical profiles exist in a given area or when the depth extent of geological units cannot be adequately constrained. Due to spatial variations of tectonic history over any given geological province, extrapolations of thicknesses over the whole area is necessarily faced with large uncertainties. Generally, interpretations will rely on rather simple models with a few layers of different thermal characteristics. Fortunately, the few kilometres of uncertainty in layer positions do not have much of an impact on geotherm calculations (Chapman, 1986).
2.2.5 Mantle heat flow
Certain additional constraints govern the depth extrapolation of heat generation values. Surface heat flow represents the sum of the heat generated in the crust plus the contribution from the mantle. Although not directly measurable, mantle heat flow is generally assumed to represent between 1/3 and 2/3 of the surface heat flow. In addition, it is natural to assume that mantle heat flow is more or less constant beneath Precambrian shields. For the Superior Province for instance, Guillou et al. (1994), assume mantle heat flow to range between 10 and 14 mW/m2. Jaupart et al., (1998) give a range of 7 to 15 mW/m2 for the Canadian Shield as a whole. The crustal contribution represents the remaining part of the surface heat flow. With these constraints it is generally possible to build plausible models of crustal heat flow.
2.3 The brittle (frictional) regime
From a rock mechanics point of view, earthquakes are dynamic slip instabilities usually occurring on preexisting sliding surfaces (Scholz, 1992). In a highly fractured intraplate environment subject to deviatoric stresses, such as the CSZ, rocks cross-cut by pre-existing fractures can slide abruptly, generating earthquakes. The process leading to fault reactivation is controlled by pressure, but is independent of rock type and temperature as a first approximation. To cause slip on a fracture (i.e. a fault in the case of an earthquake), the shear stress, induced by the existence of a tectonic stress field, has to overcome friction. Amonton's law (Jaeger and Cook, 1979) is an empirical formulation of the problem
where Tal and Sigma are the shear and normal stress acting on the plane of fracture, and µo denotes the friction coefficient. The frictional law can also be written as the Navier-Coulomb criterion:
where So is the cohesion and Pf is the pore fluid pressure, assuming compressive stresses are positive.
Assuming So = 0 and µo = 0.75, Sibson (1974) derived expressions of the critical stress difference (strength) necessary to reactivate preexisting fracture planes favourably oriented with respect to the stress field, i.e. with an angle Phi = (½) tan-1 (1/µo) approx. 26.5o with the axis of maximum compression (Sigma1)
where Sigma3 is the minimum compressive stress, Beta is a factor related to the type of faulting (3 for thrust, 1.2 for strike-slip and 0.75 for normal faults), Rho the average density of the overlying rocks at depth z, g the acceleration of gravity, and Lamda the pore fluid factor (Lamda = Pf / Rho g z).
Expanding on Sibson's formulation, Yin (1993) derived more general expressions including variations in Phi, So, and µo, as well as fracturing of isotropic rocks (with S and µ, the cohesion and friction coefficient for unfractured rocks). For thrust faulting on a preexisting fracture, the expression of the critical stress difference is
Yin (1993) also derived expressions for the limiting angles for reactivation of preexisting faults, with Phio = tan-1(µo)
If the angle between the maximum compressive axis and the fracture is outside these boundary values, a new fracture is created (if sufficient stress differences exist). Before applying these expressions to a crustal model of the CSZ, we review the main factors contained in Equations 2.9 to 2.11. Figure 2.2 shows graphically the parameters intervening in fault reactivation.
2.3.1 Crustal stress magnitude and orientation
Fault reactivation occurs when the deviatoric stresses are sufficiently large to make the shear stress overcome friction. Fault strength can be computed from equation 2.9; however, assessing the level of deviatoric stress in the crust is more difficult. In an intraplate environment such as the CSZ, there are many potential sources of deviatoric stress, some existing from plate-boundary processes and some of local origin (Hasegawa et al., 1985). In the North American plate, ridge push appears to be the main contributor to tectonic stress (Zoback, 1992a) with some contribution from basal drag. Local sources can also add to the ambient stress field. These sources are flexural stresses (such as glacial loading/unloading), lateral density contrasts, and crustal inhomogeneities (intrusions, structural boundaries, intersecting lineaments, impact craters) that can reorient and concentrate stresses. Studies of the local sources of stress have shown that glacial rebound stress are of the order of less than 5 MPa(1) (Hasegawa et al., 1985). They could act as a trigger mechanism on optimally oriented pre-existing faults on the verge of failure (Wu and Hasegawa, 1996). Lateral density variations induce stresses that do not exceed 10 MPa in the CSZ (Assameur and Mareschal, 1995). Worldwide, direct measurements and modelling of tectonic stresses suggest that the maximum crustal stress difference has an upper limit of about 100-200 MPa (see Figure 2.3 with references).
To be reactivated, faults have to be favourably oriented with respect to the stress field. With all other conditions constant, a fault with an optimum orientation will necessitate a lower critical stress difference to cause slip. In the CSZ case, small scale fractures exist at all scales due to the multiple tectonic episodes of the region. Large scale features, on the other hand, concentrate along certain trends (see section 1.3).
2.3.2 Coefficient of friction and cohesion
The coefficient of friction µ (either for intact rock µ or for fractured rock µo) is defined as (neglecting cohesion)
where Tal and Sigma are the shear and normal stress acting on the surface during sliding (Byerlee, 1978). At low pressures (less than 5 MPa), experimental data show a wide scatter as friction relates mainly to surface roughness. At higher pressures, friction becomes independent of rock types and surface roughness. Byerlee (1978) obtained two expressions for the shear stress
for Sigma3 < 114 MPa (or z < 7 km where Sigma3 is vertical), and
for 200 MPa < Sigma3 < 2000 MPa (or z > 7 km where Sigma3 is vertical).
In most cases, µ varies between 0.5 and 0.8, with a generally accepted value of 0.75. Notable exceptions are fault gouges consisting of clay minerals that can have µ as low as 0.2. Lower values for the coefficient of friction widen the range of fault orientations that can be reactivated (see Equations 2.10 and 2.11). Finally, the cohesion of the surface (or cementation strength; Sibson, 1992) is the shear stress necessary to cause slip even in the absence of normal stress. It is generally of the order of 1 to 20 MPa for preexisting fractures (Etheridge, 1983).
Recent experiments have shown that a high pressure type of brittle fracture exists where the critical stress difference increases with pressure less steeply than in the frictional case (Shimada, 1993). At higher confining pressures, the high-pressure type of fracture occurs, since compressive strength becomes lower than the frictional strength. If confirmed by further experiments, these results may provide a mechanism by which the crust's strength may not reach values as high as those predicted by Byerlee's law at mid-crustal depth. At a certain depth, the crustal strength would increase, at a lower rate than near the surface, until the ductile field is reached.
2.3.3 Pore fluid pressure
Pore fluids can favour fracturing by lowering the normal stress acting on the fracture plane and/or by chemically weakening the fault surface (Locker, 1995). The effective normal stress (Sigmaeff) is
where Pf is the pore pressure. The pore fluid factor (Lamda) can be expressed as the ratio of pore fluid pressure to lithostatic pressure
Absence of fluids implies Lamda = 0, Lamda about 0.4 (hydrostatic) means that the pores are connected from the surface to the depth considered and the pore fluid pressure is Pf = Rhow g z , where Rhow is the density of water ; and a lithostatic pressure ( Lamda about 1.0) implies that fluid pressure is equal to the overburden pressure. Suprahydrostatic pressures can occur by consolidation of sediments (in sedimentary basins) and by dehydration reactions in metamorphism (Scholz, 1990). The mechanism to sustain such pore fluid pressures, especially in areas without tectonic or volcanic activity, is still debated. There is geologic evidence that zones of suprahydrostatic pore fluid pressures are fairly common (Carter and Tsenn, 1987). As will be discussed below, this has implications for the nature of deformation at depth and for the magnitude of stress necessary to reactivate faults.
2.3.4 Density of overlying rocks
The average density of the overlying rock is used to compute the lithostatic pressure which, in the case of thrust faulting, is the minimum compressive stress (Sigma3 vertical). We use an average density Rho = 2.8 g·cm-3 for crustal rocks.
2.4 The ductile regime
The increase of temperature and pressure with depth affects the behaviour of rocks and fracture planes subject to deviatoric stress. When temperature and pressure are sufficiently high, minerals tend to flow instead of behaving in a brittle fashion. With increasing temperature, the whole rock can flow, leading to bulk flow (ductility or plasticity). Rock ductility can be modeled from extrapolations of experiments conducted at higher strain rates. The strength of the ductile part of the lithosphere is generally determined from the relation
where Sigma1 - Sigma3 is the creep strength (stress difference necessary to achieve a given steady-state strain rate Epsilon dot, R is the gas constant, T the absolute temperature and A, n and Q are material parameters which are approximately stress- and temperature-independent in the range of interest (Kirby, 1983; Ranalli, 1995).
Although strain is generally localized in highly deformed areas such as shear zones of interplate regions, it is acceptable to use average strain rate in intraplate areas where deformation affects large areas. Whereas rates can be as high as 10-12 s-1 in regions such as the San Andreas fault, in an intraplate environment such as the CSZ, Epsilon dot is assumed to be in the 10-15 to 10-16 s-1 range (Hasegawa, 1986). For comparison, 10-16 s-1 corresponds to a yearly change of 3 mm over a horizontal distance of 1000 km. There are no direct strain rate measurements for the CSZ.
The material parameters A, n and Q are known for a selection of silicate materials (examples in Table 2.1; Carter and Tsenn, 1987; Kirby and Kronenberg, 1987). The experimental uncertainties of these creep materials correspond to a 10 to 20% uncertainty in the accuracy of the creep strength (Ranalli, 1991). The creep strength of silicate rocks is approximately inversely proportional to the silica content. In the lower crust, this translates into basic rocks being stronger than acidic to intermediate rocks. Water, even in trace amounts, can enhance rock ductility (Carter and Tsenn, 1987). Thus, a number of rocks have parameters measured as dry and wet. Because finer materials are generally weaker (since creep is effected by vacancy diffusion instead of dislocation motion), small grain sizes lower rock strength and can lead to linear rheology (Ranalli, 1991). If found in bands of regional extent, these finer materials can become shear zones by concentrating strain (Kirby and Kronenberg, 1987). Section 2.6 describes how the passage from the brittle to the ductile regime can be related to the earthquake activity.
2.5 Velocity weakening-velocity strengthening transition
When the crustal temperature reaches about 300oC, the frictional properties of fault zones change and a temperature-controlled sliding stability transition occurs (see e.g. Scholz, 1990). Unlike bulk flow which is related to whole rock properties, this transition impacts on the sliding properties of the fault surface. It was found experimentally that rock volumes subjected to sliding along a fracture plane, will react differently to a change of sliding velocity, depending on the temperature. At temperatures below 300oC, rocks
Table 2.1 Creep parameters of rocks used in this study
Material | A (MPa-ns-1) | n | Q(kJ mol-1) | Reference |
Adirondack Granulite | 8.00e-03 | 3.1 | 243 | Wilks and Carter (1990) |
Orthopyroxenite, wet | 6.30e-03 | 2.8±0.5 | 271±18 | Ross and Nielsen (1978) |
Pikwitonei Granulite | 1.40e+04 | 4.2 | 445 | Wilks and Carter (1990) |
Websterite, dry | 4.00e-07 | 4.3±0.6 | 326±19 | Ave Lallemant (1978) |
Websterite, wet | 2.00e-05 | 5.3±0.4 | 382±27 | Ave Lallemant (1978) |
Clinopyroxenite, wet | 1.50e+05 | 3.3±0.9 | 490±159 | Boland and Tullis (1986) |
Clinopyroxenite, dry | 1.00e-05 | 5.3±1.1 | 380±30 | Kirby and Kronenberg (1984) |
will be in the velocity weakening mode, i.e., a change in the sliding velocity can lead to an unstable situation brought about by a decrease in the frictional shear strength. At higher temperature, velocity changes can lead to increased frictional shear strength, i.e. to velocity strengthening. The change from velocity weakening to velocity strengthening is possibly related to materials of the fault surface. For a quartzo-feldspatic crust, possible critical temperatures are 300oC and 450oC which correspond to the onsets of quartz and feldspar plasticity respectively (Scholz, 1988). Below the critical temperature, the fault is in the seismogenic velocity weakening field, while above it the fault is in the aseismic velocity strengthening field (Tse and Rice, 1986; Scholz, 1990). For granite, experiments with high fluid pressures show that the transition temperature is in the 350oC range (Blanpied et al., 1991).
2.6 Thermal and rheological controls on CSZ earthquakes
In most intraplate seismic areas of the world, the thickness of the seismogenic crustal layer correlates with surface heat flow (Sibson, 1982; Chen and Molnar, 1983; Chen, 1988). Geotherms control the maximum depth of crustal earthquakes (DMax) either by a transition from brittleness to ductile flow (the brittle-ductile transition) or by a transition from velocity-weakening to velocity-strengthening fault behaviour. Both processes occur approximately at the same temperature.
In the first model, DMax corresponds to the rheological change from a seismogenic brittle upper layer, modeled by the Navier-Coulomb criterion (Byerlee's law), to a ductile lower layer where strain is accommodated by aseismic creep (the brittle-ductile transition; see e.g. Meissner and Strehlau, 1982; Sibson, 1982; Fadaie and Ranalli, 1990). The strength of the ductile part is modelled by an extrapolation to low strain rates of steady-state creep equations (Equation 2.17). This approach has limitations (cf. Rutter and Brodie, 1992), chief among which (in addition to uncertainties in thermal and rheological parameters) is the assumption that strain is uniform, while observation shows that strain is highly localized in deep crustal exposures. However, it provides a first-order estimate of the rheology of the crust with a general rule that the critical temperature for bulk flow is about one-half of the melting temperature.
In the second model, the lower limit of seismicity is caused by a modification of slip behaviour. DMax corresponds to a temperature-controlled sliding stability transition rather than to a transition in bulk rheology (see e.g. Scholz, 1990). For a quartzo-feldspatic crust, seismicity stops when a critical temperature of 300oC for quartz-rich rocks or 450oC for feldspar-rich rocks is reached (Scholz, 1988).
In most seismic regions, both the brittle-ductile and the velocity weakening-strengthening models fit as a first approximation (which implies that the difference between the two critical temperatures is not large). However, in some areas, such as northern Switzerland (Deichmann and Rybach, 1989) deep crustal earthquakes occur, apparently in the ductile zone.
The following sections examine the rheological controls on the Charlevoix earthquake depth distribution. First, our one-dimensional geothermal modelling is described. The computed geotherms are used to infer the transition depth to bulk flow or to velocity strengthening. Finally, the frictional constraints imposed on seismicity in the mid- to lower crust are discussed.
2.6.1 Calculation of geotherms
Temperature is a critical factor in rheological models, but the uncertainties of geotherm calculations are often underestimated. Generally, these calculations use measured values of surface heat flow, and assumed values of thermal conductivity and heat production. This approach underestimates geotherm uncertainties. In the present analysis, a set of 22,000 geotherms is computed using wide ranges of values for all geothermal parameters. The calculation is based on one-dimensional heat conduction (equations 2.3 and 2.4). Fixed ranges of heat production and thermal conductivity within layers of known thicknesses are assumed.
The distribution of heat flow measurements is poor in the Grenville, with a mean value typical of Precambrian Shield regions of 41 ± 10 mW/m2 (mean value ± 1 standard deviation; Pinet et al., 1991; see Figure 2.4). There are only three surface heat flow determinations for the CSZ, one in the Appalachians (55 mW/m2) and two in the Precambrian basement: 65 mW/m2 (La Malbaie) and 57 mW/m2 (Les Eboulements; Guillou-Frottier et al., 1995). The higher than average values for the two CSZ sites are due to higher heat production rates (Guillou-Frottier et al., 1995). Two series of calculations are made in the following: one with the 41 ± 10 mW/m2 mean value for the surface heat flow in the Grenville, and one with a surface heat flow of 60 mW/m2 to consider the higher values measured in the CSZ (to be discussed later). The crust in the Charlevoix region is 43 km thick with a discontinuity at 18 km depth (from seismic refraction; Lyons et al., 1980). On the basis of geological evidence (Emslie and Hunt, 1990), the upper layer is assumed to be more felsic than the lower one which is taken to be a basic residue of partial melting, rich in plagioclase and without quartz. This hypothesis is also supported by the rather high P-velocities (6.7-6.8 km/s) measured below 18 km (Lyons et al., 1980; see Figure 5.23). In this two-layer crustal model, upper crustal heat production has been assumed to vary between 0.3 and 2.0 µW/m3 to include lateral heterogeneities such as granitic intrusives (high heat production) and large anorthositic bodies (low heat production). Lower crustal heat production has been taken to vary between 0.1 and 0.5 µW/m3 (granulite facies rocks; Ashwal et al., 1987). Thermal conductivity is assumed to vary between 2.3 and 2.7 Wm oC, a range that covers most crustal rock types. Table 2.2 lists the ranges of values for thermal parameters. Three additional conditions are used: heat production in the upper crust has to be equal to or larger than that in the lower crust, as shown by measurements of representative rocks (Ashwal et al., 1987); mantle heat flow has to be in the 10-25 mW/m2 range, which is plausible for the Precambrian Shield (it includes the 10 and 14 mW/m2 assumed for the Superior Province by Guillou et al., 1994); and thermal conductivity in the upper crustal layer is higher than or equal to that of the lower crustal layer. Computations are performed assuming all possible combinations of surface heat flow values between 31 and 51 mW/m2 (in 1 mW/m2 increments), heat production (in 0.1 µW/m3 increments) and thermal conductivity (in 0.1 W/moC increments) for the ranges specified above.
The results of 22,000 geotherm calculations can be summarized by five geotherms;
Table 2.2: Thickness and thermal parameters of crustal layers in the Charlevoix region
Thickness
(km) |
Minimum | Average | Maximum | ||
Surface | Surface heat flow (mW/m2) | 31. | 41. | ||
Upper crust | Heat production ( µW/m3 )
Thermal conductivity (W/moC) |
0.3
2.3 |
0.8
2.7 |
3.1 | |
Lower crust | Heat production (µW/m3 )
Thermal conductivity (W/moC) |
2.3 |
0.3 2.5 |
||
Mantle | Mantle heat flow (mW/m2) |
Minimum; Minimum 90%; Average; Maximum 90%; and Maximum (Figure 2.5A). Minimum and Maximum are the absolute minimum and maximum values obtained for any combination of parameters. Average represents the average of all temperature calculations. MIN90% and MAX90% are the 5 and 95 percentiles, such that 90% of the iterations fall within them. At 25 km, 90% of all computed temperatures fall between 215 and 355oC (curves MIN90% and MAX90% of Figure 2.5A). We refer to this 90% range of values when discussing seismicity.
As already pointed out by Chapman (1986), surface heat flow is the most important uncertainty in thermal calculations. In Figure 2.5B, Qmin and Qmax denote the geotherms obtained by using the lowest and the highest values of surface heat flow and the average values of the other parameters (Table 2.2). Temperature variations due to other thermal factors are constrained within an approximate ± 50oC band (arrows in Figure 2.5C). In the upper crust, variations are more influenced by heat production. While the range computed in this work is very wide (due to the inclusion of all plausible values of thermal parameters) and includes all previous values in the Precambrian Shield, the average temperature at 20 km (T 280oC) is very close to that estimated by Mareschal (1991) for the Charlevoix area. The Qmin-Qmax range also includes the most relevant geotherms of Chapman (1986) (Figure 2.5C). One-dimensional geotherm calculations that consider plausible ranges of thermal parameters clearly show the uncertainties of crustal temperatures at depth. To consider the possibility that the two surface heat flow measurements of Charlevoix are representative of a locally higher thermal regime, geotherm calculations were made with a 60 mW/m2 surface heat flow. As expected, higher crustal temperatures are obtained (Figure 2.6).
2.6.2 Ductile regime
Using equation 2.17 and the geotherms computed above, the creep strength can be calculated for given crustal material parameters. The lower crustal rocks used in this thesis are: Adirondack (felsic) and Pikwitonei (basic) granulites; dry and wet orthopyroxenites; dry and wet clinopyroxenites; and dry and wet websterites (rheological parametres given in Table 2.1). The Adirondack granulite has a much higher quartz content (21%) than plausible for the Charlevoix mid to lower crust and is only used as an upper felsic (soft) limit.
Results for three different geotherms (Minimum 90%; Average; Maximum 90%) are shown in Figures 2.7A / B / C. Geotherms lower than or equal to the Average produce a brittle crust down to 30 km or more. The maximum depth of earthquakes cannot be controlled by the brittle-ductile transition if the average Grenville geotherm applies to the Charlevoix region and the lower crust is basic. With the Maximum 90% geotherm, bulk flow may start at about 25 km, but only for the quartz-rich Adirondack Granulite. Lower strain rates could also raise the transition but the change is only a few kilometres for a decrease of two orders of magnitude. Passage to ductility at 25 km depth requires geotherms very close to the upper limit, or felsic rock composition (Adirondack granulite facies types). Ductility would be more probable with a wet lower crust as proposed by Hyndman and Shearer (1989). If the lower crust is dry as indicated by some interpretations of lower crustal petrology (Frost and Bucher, 1994), the crust would be entirely within the brittle field.
Using the geotherm computed with a surface heat flow of 60 mW/m2, bulk flow was calculated for the same materials (Figure 2.8). As expected, the brittle-ductile transition is shallower for this warmer geotherm.
2.6.3 Change from velocity weakening to velocity strengthening
In the frictional model, the transition from velocity weakening to velocity strengthening behaviour marks the lower boundary of the seismic layer. Scholz's (1988) seismo-rheological model assumes that DMax coincides with the onset of ductility of the most ductile rock component. Experimentally, 300oC marks the onset of quartz plasticity and 450oC corresponds to the onset of feldspar plasticity. Since the Charlevoix middle and lower crustal rocks are assumed not to contain quartz, the critical temperature should be the latter. The 450oC isotherm is found at about 33 km for the Maximum 90% and deeper for colder geotherms. A possible hydrolytic weakening of the feldspars could bring the critical transition temperature down to 350oC (Tullis and Yund, 1980). At 25 km, this temperature is reached by geotherms close to Maximum90% (Figure 2.4). Consequently, seismogenic control by onset of hydrated feldspar plasticity is possible with our warmest geotherms.
2.6.4 Maximum plausible stress difference
Another problem arises with the predicted stress differences necessary to cause frictional failure. The Coulomb-Navier failure criterion (with low pore fluid pressures and a static coefficient of friction; Equation 2.8 and Figure 2.6) implies very high stress differences at mid- to lower crustal depths under dry conditions. At 25 km, for example, a hydrostatic pore fluid pressure (Lamda = 0.4) and a coefficient of friction of 0.75 imply a critical stress difference for sliding of about 1200 MPa in an ideally-oriented reactivated thrust fault. These high stress differences are about one order of magnitude larger than the upper limit of 100-200 MPa discussed in Section 2.3.1.
Stress differences causing thrust fault reactivation can be computed as a function of depth, material parameters and fault orientation (Ranalli and Yin, 1990; Yin and Ranalli, 1992). The two parameters that can maintain (Sigma1 - Sigma3) below 200 MPa in the mid- to lower crust are a high pore fluid pressure and/or a low coefficient of friction. Figure 2.9A shows the critical stress difference necessary to cause thrust faulting at a depth of 20 km, assuming a static friction coefficient of 0.75, as a function of the angle Theta that the fault makes with the maximum compression direction and the pore fluid/overburden pressure ratio (Lamda). High pore-fluid pressures equal at least to 80% of overburden pressure must exist for (Sigma1 - Sigma3 ) to be less than 200 MPa. However, the steeply-dipping rift faults (70o) are not favourably oriented for thrust reactivation, assuming the sub-horizontal maximum stress direction that prevails in eastern North America (Zoback and Zoback, 1991), which may be slightly reoriented in the CSZ due to postglacial rebound (Wu, 1998). At a depth of 20 km, the limiting angle between pre-existing fault and maximum compression direction for reactivation as a thrust fault is 54o for any pore-fluid pressure below or equal to the lithostatic. An additional possibility is a reduced coefficient of friction. Figure 2.9B shows the maximum dip angle of a fault with respect to the maximum compression direction assuming that the coefficient of friction is 0.3. Reactivation of faults dipping up to 72o is now possible, and the required pore-fluid pressures are considerably less. Assuming a 200 MPa maximum stress difference, it is possible to map the relationship between the pore-fluid pressure, the coefficient of friction and the angle between Sigma1 and the fault surface (Figure 2.10). There is also the possibility that the paleo-rift faults become listric at depth. According to Figure 2.10B, fractures that make an angle of 40o with the axis of maximum compression can be reactivated with coefficients of friction in the 0.6-0.75 range. Although high pore-fluid pressures are still required (Lamda of about 0.8), listric faults could be reactivated with standard values for the coefficient of friction. Finally, the high pressure type of fracture may also account for the maximum crustal strength at a level slightly above 200 MPa (Shimada, 1993; Figure 2.11).
2.7 Discussion and conclusions
In Charlevoix, D99%, D80% and DDec do not appear to correspond to a transition to ductile bulk rock behaviour unless the geotherm is near the maximum possible value and/or the lower crust is quartz-rich. This conclusion is based on one-dimensional geothermal modelling of Grenville surface heat flow for a wide range of thermal parameters and on an examination of the rheology of materials likely to form the lower crust. None of the materials can be ductile within our 90% confidence limits for the geotherms. Only a high surface heat flow ( 60 mW/m2) can lead to a passage to ductility at D99%, with a quartz-rich mid- to lower crustal composition. However, the high heat flow value (Q 65 mW/m2 ; Guillou-Frottier et al., 1995) appears to be affected by the high thermal conductivity (J.-C. Mareschal, pers. comm.).
Therefore, although there is a global correlation between surface heat flow and maximum depth of crustal earthquakes (Sibson, 1982; Meissner and Strehlau, 1982; Chen and Molnar, 1983), a detailed analysis of the CSZ shows that ithe brittle-ductile transition is probably deeper than the D99% of earthquakes. It should be noted that the present model uses a basic lower crust while most models assume a uniform felsic crust with the rheology of quartzite or granite. The present choice, on the other hand, is based on inferences on the composition of the mid- to lower crust of the Charlevoix region. Another difference between this approach and models of other regions is the consideration of temperature uncertainties. Most geotherm models assume fixed values for thermal parameters.
The passage from velocity weakening to velocity strengthening frictional behaviour may explain D99% . Feldspar only starts flowing at about 450oC, but this temperature can be lowered to about 350-400oC by the presence of H2O (Tullis and Yund, 1980). Consequently, at 25 km depth, hydrated feldspar can possibly become ductile with the warmest geotherms (including the one with 60mW/m2 surface heat flow). Below D99% , the absence of three local instability factors may explain the quasi-complete absence of earthquakes. The available stress difference may not reach the frictional strength level. High pore-fluid pressures may exist only above D99% , possibly due to hydration reactions that eliminate any free water at temperatures over 250oC (Frost and Bucher, 1994). Finally, below D99% fractures may be less numerous or more stable possibly due to healing. Figure 2.12 presents a synoptic model of seismogenic factors in this region.
Unless the high-pressure type of fracturing is validated by further testing, high pore fluid pressure and/or low coefficient of friction must exist not to exceed the maximum crustal stress difference level. Earthquakes occur where these conditions exist. The presence of deep crustal fluids and of high pore fluid pressures has gained much support in recent years (Hyndman and Shearer, 1989; Meissner and Wever, 1992). Examples for eastern Canada are numerous: Zoback (1992b) used a high pore-fluid pressure to explain the apparent mismatch between the regional compressive stress direction and the focal mechanism of the 1979 mbLg 5.0 Charlevoix earthquake; shear wave splitting observed for Charlevoix earthquakes can be explained by saturated cracks (Buchbinder, 1989); in the Miramichi earthquake area in New Brunswick, near-lithostatic pore-fluid pressures can explain reactivation of high angle faults (Sibson, 1989). One open question, however, is how such high pore fluid pressures are created and maintained. In the CSZ, no active geological process, such as tectonism, sediment compaction or metamorphic dehydration, are likely to give rise to suprahydrostatic fluid pressures. In addition, it has been argued that at temperatures higher than 250oC, free water cannot exist due to a rapid consumption by hydration reactions (Frost and Bucher, 1994). Although the debate about the presence or absence of mid- to lower-crustal fluids is open, high pore-fluid pressures appear to solve the seismological problem of maximum possible stress difference at mid- and lower-crustal depth.
The alternative is a low friction coefficient. In a dry thrust environment with sub-horizontal maximum compressive stress, a low friction coefficient is necessary to reactivate steeply-dipping faults. This possibility was also suggested by Zoback (1992b). Near the surface, low friction coefficients may correspond to thick fault gouges made up of clay minerals; at depth, however, these minerals lose their lubricating properties due to dehydration (Scholz, 1990). A low friction coefficient at mid-crustal depth may correspond to highly fractured zones. Variations in hypocentral depths along the trend of the St. Lawrence River, and variability in reactivated fault planes shown by focal mechanisms (Figure 1.24) indicate that Charlevoix earthquakes occur in zones of intense fracturing.
In terms of the seismic hazard estimate in the St. Lawrence valley, the presence of paleo-rift faults alone does not explain earthquake occurrence. Local fault conditions, such as transient enhanced pore-fluid pressures, leading to frictional instabilities should also be considered.
Worldwide, large intracontinental earthquakes tend to occur at the base of the seismogenic layer defined by micro-earthquakes, which represents the peak in the shear resistance of a rheological profile (Sibson, 1984). In contrast, the two largest CSZ earthquakes for which focal depths are known (1925 mb 6.2: Bent, 1992; and 1979 mbLg 5.0: Hasegawa and Wetmiller, 1980) occurred at about 10 km depth. In addition, the 5 to 15 km depth range also includes magnitude 3.0 earthquakes for the period 1977 to 1997. Based on our rheological modelling, the concentration of larger earthquakes between 5 and 15 km depth cannot be associated with the much deeper brittle-ductile transition. Again, local factors, such as high pore fluid pressures, degree of fracturing, or lower coefficient of friction, may explain the enhanced activity of this crustal level.
The mbLg 6.5 1988 Saguenay earthquake occurred at 29 km depth, in the Grenville Province. Based on our rheological modelling, it probably occurred in the brittle crust, but much closer to the brittle-ductile transition, if a warm Grenville geotherm (Max90%) is considered (a situation similar to Figure 2.8). In this case, the earthquake may have occurred in the strongest part of the crust, as proposed by Sibson (1984). This would be consistent with the unusually high stress drop of this event (Atkinson and Boore, 1998).
If the two higher than average CSZ surface heat flow measurements of Guillou-Frottier et al. (1995) are representative of a regional anomaly in the Grenville, the D99% may correspond to the brittle-ductile transition (assuming the warmest geotherm). In this case, many earthquake concentrations in the 20 to 30 km depth range could be seen as "sticking patches" in a surrounding ductile regime (e.g. Sibson, 1984). Interestingly, most earthquakes below 20 km depth occur in an elongated zone, that could be located between more ductile areas. It would be interesting to see if the next large CSZ earthquake nucleates at 20 to 30 km depth, near one of these locking points.
Although this study examines the case of the Charlevoix seismic zone, most questions raised in this analysis also apply to other seismically active areas of the continental crust. It appears that the onset of bulk flow is only one factor affecting the maximum depth of intraplate earthquakes. Charlevoix earthquakes are not controlled by the brittle-ductile transition. The probable composition of the mid- to lower crust and the surface heat flow show that rocks are still brittle at D99%. A passage from velocity-weakening to velocity-strengthening at 300-350oC may correspond to D99%, but only if this transition temperature applies to hydrated basic rocks, where feldspar is the most ductile material. Some instability factors may explain the occurrence of earthquakes: high pore-fluid pressure and low friction coefficient, possibly related to unhealed zones of intense fracturing. Independently of the pore-fluid pressure, a low coefficient of friction may be another factor to account for the thrust reactivation of steeply-dipping faults.
Figure 2.1 Focal depth distribution of the Charlevoix earthquakes for the period October 1977 to December 1997. Events are grouped in one kilometre intervals. Two curves are drawn: one for absolute numbers (in light blue), and one for the cumulative percentage (in grey). See text for description of DDec ; D80% ; D99% , and Dmax .
Figure 2.2 Sketch representing the various factors intervening in fault reactivation.
Figure 2.3 Distribution of inferred maximum stress differences. The numbers refer to: 1) Ridge push: Zoback and Magee (1991); 2) Upper crustal amplification of ridge push: Kusznir and Bott (1977); Hasegawa et al. (1985); 3) Tensile fracturing : Etheridge (1983); 4) Earthquake stress drops in intraplate environments: Hasegawa et al. (1985); 5) Estimate along the San Andreas fault: Lachenbruch and Sass (1980); 6) Hydrofracturing-strain relief: McGarr and Gay (1978); 7) Composite estimates from the summation of eastern Canadian stress sources: Hasegawa et al. (1985); 8) Isostatic considerations: Molnar and Lyon-Caen (1988); 9) Mylonitization: Christie and Ord (1980); 10) Paleopiezometers: Kohlstedt and Weathers, 1980); 11) KTB borehole value measured at 6 km: Zoback et al. (1993); 12) Estimates along the Alpine fault of New Zealand: Scholz et al. (1979); 13) Mechanical twinning: Tullis (1980). Except for mechanical twinning, a 200 MPa is a conservative upper bound estimate of the maximum crustal stress difference.
Figure 2.4 Distribution of heat flow measurements in the Grenville Province, with circles proportional to the surface heat flow values, given in mW/m2 (from Jessop et al., 1984; Drury 1987; Pinet et al., 1991; Guillou-Frottier et al., 1995).
Figure 2.5 Computed geotherms for the Grenville region, (A) using parameters of Table 2.2. Five curves are shown: the absolute Maximum (MAX) and Minimum (MIN) geotherms; the Average (Av) geotherm which represent the average value of all geotherms at a given depth; and the MIN90% and MAX90% which enclose 90% of all geotherms. Ninety-nine percent of all earthquakes occur above D99% . (B) Temperature variations due to one parameter while others are kept at their average values of Table 2.1. Surface heat flow is the most important factor, while all others constrain geotherm fluctuations within 50oC when the average heat flow is used. (C) Comparisons of MIN90%, AV, and MAX90% geotherms with the geotherms of Chapman (1986) for a surface heat flow of 40 mW/m2 (Ch40) and 50 mW/m2 (Ch50). Both fall within the present 90% confidence margin.
Figure 2.6 Geotherms calculated with a surface heat flow of 60 mW/m2 (Min60: minimum; Ave60: average; Max60: maximum), compared with the minimum 90% and maximum 90% geotherms of Figure 2.5A.
Figure 2.7 Crustal strength as a function of depth. The brittle part is computed from the Ranalli and Yin (1990) formulation using µ = 0.75, S0= 10 MPa, Rho= 2.6 g·cm-3 and Lamda = 0.0 (dry) and 0.9. The ductile part is computed using equation 2.17, with rheological parameters given in the text, a strain rate of 10-16 s-1, and three geotherms: A) Minimum 90% ; B) Average; C) Maximum 90%. The rocks used are: Adirondack Granulites (AG); wet Orthopyroxenite (wO); Pikwitonei Granulite (PG); wet Websterite (wW); and wet Clinopyroxenite (wC). Rocks that do not appear in the ductile part of the diagrams have a creep strength higher than 10000 MPa. MaxS placed at 200 MPa indicates the maximum crustal stress difference discussed in the text and interpreted from Figure 2.3.
Figure 2.8 Crustal strength as a function of depth, with same creep parameters as Figure 2.7 except for the geotherm which is the Max60 geotherm of Figure 2.6.
Figure 2.9 Stress difference (MPa) necessary to reactivate a fault located at a depth of 20 km for various angles between the maximum compressive stress axis and the fault and various ratios of pore fluid pressure Lamda: A) With a friction coefficient µ = 0.75 (standard); B) with a low friction coefficient µ = 0.3. In A), a maximum stress difference of 200 MPa can only be reached with very high Lamda. With a low friction coefficient, the maximum reactivation angle is increased and the necessary pore-fluid pressure is decreased.
Figure 2.10 Map of the relation between pore-fluid pressure, coefficient of friction and reactivation angle, assuming a maximum stress difference of 200 MPa. A) Three dimensional view, with colours representing the pore fluid pressure ratio; B) a two dimensional view. Taking an angle of 30o as an example, one can see that the assumed value of 0.75 for the friction coefficient requires pore fluid pressures in the > 0.80 range. Similarly, if the angle between Sigma1 and the fault plane is 70o, a < 0.35 friction coefficient is necessary to reactivate the fault. Areas in black are outside the range 0.0 <= lamda <= 1.0.
Figure 2.11 Example of the impact of a change in fracture type for the rheological profile. At a given depth (10 km in this example), the brittle fracture changes from Byerlee's law (A), between the dry (black) and the wet (blue) curves, to a high pressure type of fracturing (A'). This change in fracture type maintains the strength near the maximum plausible crustal level (MaxS). Earthquakes can occur down to the ductile level without the need for very high stress difference.
Figure 2.12 Synoptic model of the factors controlling the earthquake occurrences in the CSZ. On a cross-section perpendicular to the St. Lawrence River, the earthquake distribution is shown, with D99% at 25 km depth. The major factor influencing the depth distribution is the temperature. Two possibilities exist: an average geotherm, where the 300oC is reached at about 28 km, or the Max90% geotherm, where the 300oC is reached at about 20 km. D99% may correlate with a number of factors: the passage to velocity strengthening fault behaviour; a lessening of the pore-fluid pressure; or an increase in the friction coefficient. The brittle-ductile transition however lies below D99% . In the brittle mid-crust, earthquakes occur in regions of enhanced pore-fluid pressures or of low coefficient of friction, possibly due to highly fractured zones.
1. One MPa corresponds to 10 bars.