## Abstract

Plane strain thermo-mechanical finite-element model experiments are used to investigate the effects of frictional–plastic strain softening and inherited weakness on the style of lithospheric extension. The model results are compared with the Newfoundland–Iberia conjugate rifted margins with the goal of understanding the lithospheric properties that controlled their evolution during rifting. Our proposition is that coupling between the plastic–viscous layering, acting together with frictional–plastic strain softening localized on inherited weak heterogeneities, can explain the initial wide rift and distributed rift basins that are later abandoned in favour of a narrow rift in which mantle lithosphere is exhumed to the surface. The models comprise uniform composition viscous and plastic layers in which focused deformation is nucleated on either a single weak ‘seed’ or a statistical white noise distribution of inherited strain. Strain softening of frictional–plastic layers acts as a positive feedback mechanism that creates localized shear zones from the inherited weak heterogeneities. The sensitivity of deformation to the choice of softening parameters and the type of inherited noise is examined in cases where the deeper part of the crust is either weak or strong.

Lithosphere-scale models with a single weak seed exhibit a range of asymmetric and symmetric rifting modes that are mostly determined by the feedback between two primary controls, coupling between the plastic and viscous layers and strain softening. Decreasing and increasing the rifting velocity can change the mode, and asymmetry is strongest in models with low rifting velocities and a strong lower crust. Analysis of equivalent simple-bonded plastic–viscous two-layer models using the minimum rate of dissipation principle demonstrates that the mode selected depends on the division of the dissipation between the layers. Criteria developed on minimizing the total dissipation show how mode selection changes with increasing viscosity, or rifting velocity, from the: asymmetric plug or half-graben (AP) mode; through the symmetric plug or graben (PS) mode, to the distributed pure shear (PS) mode. Numerical models confirm these results.

Models with statistical white-noise-inherited strain have similar modes to those with a single seed. In addition, modes with multiple sets of shear zones develop in the plastic layer for a range of intermediate parameter combinations. We believe that distributed noise in combination with a weak lower crust and slow extension can produce model results in accord with general features of the Newfoundland–Iberia conjugate margins; an initially distributed wide rift mode, followed by a late-stage narrow rift with a significant component of mantle exhumation.

Despite the large number of studies of continental rifts and rifted continental margins, many of which have produced high-quality data (Keen *et al.* 1987*a*, *b*, 1989; Mutter *et al.* 1989; Boillot *et al.* 1992; Sibuet 1992; Brun & Beslier 1996; Louden & Chian 1993; Dean *et al.* 2000; Funck *et al.* 2003, 2004; Hopper *et al.* 2004), we still lack a unified understanding of the mechanical and thermal processes that control their extensional geometry. This problem is compounded by the wide range of styles that include both non-volcanic and volcanic types, in which the latter exhibits direct evidence of voluminous mafic extrusives, and/or indirect evidence of intruded and underplated mafic magmas. In addition, the relative importance of active v. passive rifting, that is whether rifting is driven by active mantle upwelling or whether mantle upwelling is a passive response to lithospheric extension (Salveson 1978; Sengor & Burke 1978; Turcotte & Emmerman 1983; Ruppel 1995; Huismans *et al.* 2001) remains to be determined.

Important contributions to the problem of passive rifting leading to non-volcanic rifted margins defined possible end-member styles of lithospheric extension and their associated geometries, pure shear (McKenzie 1978), simple shear (Wernicke Burchfield 1982; Wernicke 1985) and combinations of these styles, for example Lister *et al.* (1986), referred to here as compound models (Fig. 1a–c). These kinematic models subsequently served as templates for the interpretation of observations. A primary focus in these interpretations has been whether rifting is symmetric, corresponding to the pure shear template, or asymmetric as predicted by the simple shear template (Fig. 1a, b).

More generally, our lack of understanding of passive rifting at non-volcanic continental rifts and rifted margins can be expressed in terms of the following questions.

Can the geometrical styles of natural continental rifts and rifted margins be classified according to modes, such as the narrow, wide and core complex rift modes proposed by Buck (1991) (Fig. 1 a–d).

If so, how do the kinematic template styles of extension listed above, and their associated symmetries and asymmetries (Fig. 1), relate to the modes proposed by Buck (1991) (Fig. 1)?

If natural rifts can indeed be classified according to modes of extension, what controls the selection of these modes?

Conversely, if natural rifts exhibit various geometries that defy any simple mode classification, what controls this variability?

Numerical modelling has the potential to provide at least partial answers to these questions. Model experiments can be designed to test the role of a range of factors and processes on the styles of lithospheric extension. If the models are reasonable, albeit simplified, representations of their natural counterparts, the model results can indicate which processes are most likely to be important and how processes compete and tradeoff in different circumstances.

We have chosen to focus on aspects of the rifting problem that we believe were important in the development of the Newfoundland–Iberia conjugate margin system and we address this specific application at the end of the paper. More generally, observations confirm that faulting and deformation on shear zones are important during rifting. At the small scale, deformation on faults and shear zones is inherently asymmetric. However, the question of the role of shear zones on the large-scale asymmetry of passive margins is still unanswered. This is partly a result of the lack of good constraints on the structure of the continental mantle lithosphere and its evolution during extension. Nevertheless, data do point to asymmetric geometries of conjugate margins (Boillot *et al.* 1992; Sibuet 1992; Louden & Chian 1999), and, for example, large-scale low-angle detachments and features interpreted to be high-angle shear zones that extend into the upper mantle (Brun & Gutscher 1992; Whitmarsh *et al.* 2001). There is some indication of patterns, or modes, in margin geometries. The significant variability in crustal-scale geometry that commonly occurs over quite short distances along strike on many margins poses an additional problem (see, for instance, the Iberia–Newfoundland margin system: Keen *et al.* 1987*a*, *b*; Torne *et al.* 1994; Funck *et al.* 2003; Dean *et al.* 2000; Pérez-Gussinyé *et al.* 2003). Either this strike variability, or margin segmentation, is an inherent and perhaps predictable consequence of the mechanics of rifting in three dimensions, or it represents an effect of inherited inhomogeneous lithospheric properties, ‘noise’, in the system. For these reasons our general focus in the models described here is on processes that create shear zones and lead to mode selection, and on the effects of noise.

We first assume that the lithosphere can be represented by a prototype uniform-layered model comprising frictional–plastic and viscous layers. We then consider the effects of different factors or processes on the style of extension of this model lithosphere. This approach is designed to answer the following particular questions.

What is the effect of strain-softening of the frictional–plastic parts of the system on the geometry during extension and does this feedback into mode selection?

Can the mode selection be explained by a general governing principle, specifically the minimum rate of dissipation of energy?

What is the effect of inherited heterogeneity, ‘noise’, in the model crustal properties on the style of extension? Does mode selection still occur in the presence of noise, or does noise lead to incoherent extension that is merely an expression of the noise?

In section on ‘Thermo-mechanical models of lithosphere extension’ it is shown that strain-softening leads to mode selection, and to asymmetric and symmetric styles of lithospheric extension. By analysing simple two-layer models we then demonstrate in the section on ‘Analysis of simple models using the principle of minimum energy dissipation’ that the minimum dissipation principle can be used as a guide to mode selection, at least in the numerical models. The section on ‘systems with statistical heterogeneity’ contains a preliminary treatment of the effects of distributed noise in the system on mode selection, and on how noise can lead to wide and narrow styles of lithospheric extension. Lastly, in ‘implications for the Iberia–Newfoundland conjugate margin system’ we ask whether what has been learned from the numerical models can be applied to the Newfoundland–Iberia conjugate margins to determine the controls on the tectonic style of this system. These sections are prefaced by brief reviews of processes that result in strain-dependent strength of rocks, the role and effect of inherited heterogeneity, and the two ways, deterministic and statistical, that this heterogeneity is introduced into the models. The section on ‘Thermo-mechanical numerical model’ describes the numerical models.

### Strain-dependent strength of rocks and its representation in the numerical models

The strength of faults and shear zones in the crust and mantle lithosphere on geological timescales is a controversial subject. They are interpreted to be either much weaker (Chéry *et al.* 2001; Provost & Houston 2003) or of similar strength (Scholz 2000) to adjacent relatively non-deforming regions. Here, we accept that faults may be very weak and outline mechanisms that produce weak faults and shear zones. We then describe a simple parametric representation of these processes that we use in the numerical models.

The frictional, or brittle, strength of dry rocks is represented to a first approximation by the Coulomb yield criterion for incompressible deformation in plane strain:
1
where (*J*′_{2})^{1/2} is the second invariant of the at yield deviatoric stress tensor, *c* is the cohesion, ϕ is the internal angle of friction and *p* is the dynamic pressure (mean stress).

Several mechanisms may contribute to strain softening, that is a reduction of the strength of the rock (deviatoric stress at yield, equation 1), with increasing strain but their relative and absolute importance is poorly constrained. At the scale of the crust, the brittle strength is dominated by the pressure-dependent term (equation 1). Consequently, cohesion loss results only in a limited reduction of the strength of the brittle layer in the order of 2–10% depending on the ambient dynamical pressure and the depth of the brittle–ductile transition in the crust (Fig. 2a). Other mechanisms, such as formation of gouge and mineral transformations which reduce ϕ with increasing deformation, can potentially result in much larger frictional strength reduction. For example, experimental work (Bos & Spiers 2002) on quartz–feldspar aggregates and their analogues has shown that mineral transformations in granitic rocks in the presence of water may produce mica-rich shear zones with very low effective internal angles of friction (ϕ=3°–10°), resulting in a strength reduction of the order of 50–80% (Fig. 2b). Similar deformation- assisted mineral transformation mechanisms have been documented in mafic lower crustal and upper mantle lithosphere rocks in shear zones in the Ivrea zone, which formed during extension (Handy & Stuenitz 2002).

Transient and long-term fluid pressure variations may be significant factors in weakening brittle rocks (Sibson 1990; Rice 1992; Ridley 1993; Streit 1997; Ingebritsen & Manning 1999; Connolly & Podladchikov 2000). Pore fluid pressure, *p*_{f}, acts to reduce the dynamic pressure and *p* in equation 1 is replaced by the effective pressure, *p*−*p*_{f}. It can be seen that the rock strength can become very small when *p*_{f}∼*p*. An approximate way to consider the effect of pore fluid pressures is to introduce an effective internal angle of friction ϕ_{eff} such that *p* · sinϕ_{eff}= (*p*–*p*_{f}) sinϕ=*p*(1–λ)sinϕ, where λ=*p*_{f}/*p* is the pore pressure ratio. (Note that this definition differs from the Hubbert–Rubey pore pressure ratio for which λ=*p*_{f}/*p*_{overburden}, where *p*_{overburden} is the weight of the overburden.) For example, using this definition of ϕ_{eff}, rocks with ϕ=30° have ϕ_{eff}≈18° in crust with hydrostatic pore fluid pressure. When the full range of pore fluid pressure variations that may occur during tectonic deformation is taken into account ϕ_{eff} may vary from dry, ϕ_{eff}≈30°, through the hydrostatic value to strongly overpressured, ϕ_{eff}≈0° (Fig. 2c). We use this effective internal angle of friction approach and adopt reference initial values of ϕ_{eff}=30°, 15° and 7° in order to investigate the sensitivity of the model results to a range of pore fluid pressure regimes.

We also include a simplified parametric form of strain-dependent weakening in the numerical model experiments by linearly decreasing the internal angle of friction with the second invariant of the deviatoric strain, (*I*′_{2})^{1/2}. In most of our models strain softening occurs in the range 0.5<(*I*′_{2})^{1/2}< 1.5. By using this approach we avoid the details of the physical mechanisms of strain softening, which are not well known in general but we can determine the sensitivity of the models to a range of possible strain-softening mechanisms.

Strain-rate dependent weakening has also been used (Behn *et al.* 2002) to represent the rate dependence of frictional strength, which is observed in laboratory experiments. This approach characterizes the frictional behaviour of faults on seismic timescales (Scholz 1990). In contrast, our modelling addresses fault strength variations on geological timescales; therefore, we focus on the strain dependence of material properties.

### Role of inheritance and heterogeneity

We believe that the style of deformation of the lithosphere during rifting reflects the interplay of extension of pristine lithosphere with the superimposed effects of inherited structural heterogeneities, which can include both anomalously weak and strong regions. Models of extension of pristine lithosphere will display deformation modes that reflect the continuum nature of the system. Uniform pure shear is one such mode. However, the observed styles of rifts and rifted margins are highly variable, suggesting that inherited heterogeneities, ‘noise’ in the system, have a strong influence on its deformation. We cannot argue that the lithosphere is pristine. Therefore, it is critical to understand whether the noise overprints the inherent modes, or whether the underlying modes persist despite the noise. The noise is even more important if it contributes to the positive feedback effect of strain weakening, as discussed in the previous section. A small inherited weak heterogeneity, which is initially the focus of mildly enhanced strain, can become the seed for major localized deformation through this feedback mechanism.

We consider two types of interaction between inherited ‘noise’ and strain softening. In sections on ‘Thermo-mechanical models of lithosphere extension’ and ‘Analysis of simple models using the principle of minimum energy dissipation’ we focus on the behaviour of deterministic model systems that include a single weak spot, or seed. This seed is taken to represent the weakest single zone of inherited damage that controls the strain localization process and, therefore, dominates over all other pre-existing heterogeneities. A sensitivity analysis enabled us to conclude that weak seeds need to have a critical size and are most effective in localizing deformation when placed in the strongest part of the system, the upper mantle lithosphere in our models, thereby causing the greatest strength contrast with the surrounding material. From this analysis, we conclude that the system will preferentially select this type of seed over others, and we therefore focus on seeds placed in the upper mantle lithosphere.

The approach that employs a single deterministic weak seed fails to consider that the weakest point may be ignored in favour of a distribution of inherited weak heterogeneities that link together to form even weaker shear zones. In the sections on ‘Systems with statistical heterogeneity’ and ‘Implication for the Iberia–Newfoundland conjugate margin system’ we consider the behaviour of systems characterized by a statistical distribution of inherited weaknesses that both compete and co-operate by linking to form shear zones. In this type of system, the continuum deformation activates the noise to create nascent-linked shear zones. Multiple shears may develop, evolve and compete. These two types of initial weakness, a single deterministic location or a statistical distribution, can be interpreted to be equivalent to natural lithospheric systems that are either dominated by one major zone of weakness or an evenly distributed background level of ubiquitous damage in the form of minor faults or shear zones that have not developed into penetrative shears at the crust or lithospheric scale.

## Thermo-mechanical numerical model

### Model description

We use an Arbitrary Lagrangian–Eulerian (ALE) finite-element method for the solution of thermo-mechanically coupled, plane-strain, incompressible viscous–plastic creeping flows (Fullsack 1995; Willett 1999; Huismans & Beaumont 2003). The model is used to investigate extension of a layered lithosphere with frictional–plastic and thermally activated power-law viscous rheologies (Fig. 3). For details on the numerical method see Fullsack (1995), Willett (1999) and Huismans & Beaumont (2003).

When the state of stress is below the frictional–plastic yield the flow is viscous. We use both linear Newtonian viscous materials (see later sections in this paper) and temperature-dependent non-linear power-law rheologies (see later sections in this paper) based on laboratory measurements on ‘wet’ quartzite (Gleason & Tullis, 1995) and ‘dry’ olivine (Karato & Wu 1993), where the latter also includes the pressure dependence of viscosity. The effective viscosity, η, in the power-law model is of the general form:
2
where (*İ*′_{2})^{1/2} is the second invariant of the deviatoric strain rate tensor (½ε̇′_{ij}ε̇′_{ij})^{1/2}, *n* is the power-law exponent, *A* is a scaling factor, *Q* is the activation energy, *V* is the activation volume (which makes the viscosity dependent on pressure, *p*), *T* is the absolute temperature and *R* is the universal gas constant. *A*, *n*, *Q* and *V* are derived from the laboratory experiments, and the parameter values are listed in Table 1.

The reference parameter values for wet quartz listed in Table 1 lead to a weak viscous lower crust. Models of this type are described as weak crust models. In other models crustal viscosity, η(wet quartz), is everywhere increased by a factor of 100 to achieve a crust that is totally in the frictional–plastic regime. Models of this type are termed strong crust models. The viscosity scaling represents a simple technique that creates either strong frictional lower crust or moderately weak viscous lower crust without recourse to additional flow laws, each with its own uncertainties. The scaling can either be interpreted as a measure of the uncertainty in the flow properties of rocks where flow is dominated by quartz (e.g. ‘wet’ or ‘dry’) or to represent strong lower crust dominated by minerals that deform in the frictional regime for the conditions (temperature, strain rate) chosen for the model experiments.

The plastic (frictional or brittle) deformation is modelled with a pressure-dependent Drucker-Prager yield criterion which, with suitable adjustment of constants, is equivalent to the Coulomb yield criterion for incompressible deformation in plane strain. Yielding occurs when:
3
where (*J*′_{2})^{1/2}=(½ σ′_{ij} σ′_{ij})^{1/2} is the second invariant of the deviatoric stress, *c* is the cohesion, ϕ_{eff} is the effective internal angle of friction, defined earlier. With appropriate choices of *c* and ϕ_{eff} this yield criterion can approximate frictional sliding in rocks. Plastic flow is incompressible. Strain softening is introduced by a linear decrease of ϕ_{eff} (ε) over a range of strain where ε represents the second strain invariant, (*I*′_{2})^{1/2} (Fig. 3, Tables 1 and 2). We have tested the dependence of the model results on mesh resolution and on the range of strain over which weakening occurs. Increasing the mesh resolution leads to narrower shear zones and to earlier localization because strain accumulates at a higher rate. Moderate changes in resolution do not, however, alter the main character of any particular model result. Also, the models do not show a strong sensitivity to the range of strain over which softening occur. The tests indicate that the major control on the model behaviour are the threshold at which strain softening starts and the amount of softening that occurs. The strong positive feedback between softening and strain accumulation implies that once initiated strain softening occurs rapidly.

In addition to solving the equilibrium equations for plastic–viscous flows, we also solve the thermal evolution in two dimensions. The mechanical and thermal systems are coupled though the temperature dependence of viscosity and density, and are solved successively during each time step. Initial conditions and other thermal properties are given in Figure 3 and Table 1.

## Thermo-mechanical models of lithosphere extension

We start with an overview of a more comprehensive set of models from Huismans & Beaumont (2002, 2003). The models selected here have frictional-plastic strain softening, and deformation is nucleated by deterministic noise in the form of a single weak seed (Fig. 3). Most of the results can be classified as narrow rifts as defined by Buck (1991). We review the sensitivity of deformation style to the extension velocity and the strength of the middle and lower crust.

### Sensitivity of model results to rifting velocity

In order to demonstrate that rifting velocity is an important control on the extension mode we show results for a reference velocity of 0.3 cm year^{−1} and for end-member variations of velocity of 0.06 and 10 cm year^{−1} (Fig. 4). The models have strong middle and lower crusts, and include strain softening of the frictional–plastic rheology.

At the reference velocity (Fig. 4b) strain softening in the frictional-plastic parts of the lithosphere results in strong localization in a single system of ‘faults’ (strictly shear zones), resulting in an asymmetric mode of extension. Middle and lower crust has been tectonically cut out and mantle lithosphere has been exhumed by large-scale frictional detachment. The proto-margins show distinct differences resembling an upper and lower plate conjugate margin pair. The main cause of asymmetry in the model is the frictional-plastic strain softening. After an initial phase of symmetric extension, feedback of strain softening results in preferential weakening of one of the two conjugate shears. The single weak shear zone remains the weakest part in the system and results in lithosphere-scale asymmetry. During later stages thermal advection produces a hotter, ductile lithosphere and a viscous necking style, and the final breakup phase is symmetric. However, asymmetry remains from the initial stages of the model evolution.

At high rifting velocities (Fig. 4a) the style of extension is markedly different. The proto-margins are essentially symmetric. We have interpreted this model behaviour in terms of the increased role of viscous rheologies in the system and a stronger viscous coupling, where higher strain rates equate to higher viscous stresses at the base of the frictional layer (Huismans & Beaumont 2002, 2003). The viscous layer promotes a distributed, symmetric style of extension, and the tendency for localization and asymmetry given by the frictional–plastic strain softening is suppressed.

The effect of a low end-member value for the extension velocity, *v*=0.06 cm year^{−1}, is illustrated in Figure 4c. The model is strongly asymmetric throughout its evolution. At this velocity thermal conduction is more efficient than thermal advection (Peclet number <1). Consequently, the thermal evolution in the model lithosphere is close to the conductive limit, and the rift zone is relatively cool and deforms in the frictional regime. The frictional weak shear zone therefore remains the single major weakness in the lithospheric system that allows for ongoing asymmetric extension.

### Sensitivity of model results to strength of the lower crust

We next examine the sensitivity of the model behaviour to the strength of the lower crust. The same frictional strain softening occurs as in the reference model, and in addition viscous flow following a wet quartz flow law now occurs in the weaker lower and middle crust (Fig. 5a). During the initial stages the crustal asymmetry is diminished because the conjugate frictional shears sole out in the weak ductile lower crust. Asymmetry in the lower lithosphere is still apparent in the shear zone cutting the mantle lithosphere. At later stages a set of synthetic frictional-plastic crustal shear zones form on either side of the central rift zone, thereby creating blocks bounded by frictional shear zones that sole out in the lower crust and that propagate beneath the rift flanks. During the late stage of deformation thermal advection results in nearly symmetric ductile necking of the lower lithosphere. The reference model result for a strong lower crust is given in Figure 5b.

### Factors controlling mode selection and asymmetry during extension

These results indicate three fundamental rift modes for the type of model described above: (1) lithospheric-scale symmetric rifting; (2) lithospheric-scale asymmetric rifting; and (3) asymmetric upper lithosphere rifting concomitant with symmetric lower lithosphere rifting. Strain softening is the fundamental cause of the asymmetry and in a purely frictional-plastic model this would always lead to a phase of asymmetrical extension. In more complex models that include a viscous lower crustal layer the degree of asymmetry accompanying strain softening also depends on the rate of extension and on the viscous strength of the lower crust.

The following factors need to be considered in regard to the efficiency of the strain-softening feedback in the frictional-plastic layer.

Geometrical hardening (Fig. 6d) may occur as the model evolves, for example by the rotation of shear zones, so that the active shears are no longer geometrically viable. This type of hardening opposes the continued increase in asymmetry and, depending on the amount and range over which strain softening occurs, even very weak shear zones will lock and be abandoned.

The extent to which the distribution of frictional-plasticity dominates the model lithosphere. A range of behaviours exists between models where by frictional-plasticity is the dominant process in the extending lithosphere, which exhibit the most asymmetry, and those where geometrically larger ductile regions promote symmetry.

The relative integrated strengths, as opposed to geometries, of the plastic and ductile regions. For example, ductile strength increases with strain rate, producing higher stresses at higher rifting velocities, whereas the frictional strength is independent of strain rate. This strain-rate dependence of strength explains the velocity sensitivity (Fig. 6d). The strong ductile region suppresses the asymmetric deformation preferred by the overlying frictional-plastic layer; ductile symmetry dominates over plastic asymmetry (Buck

*et al.*1999, 2003; Huismans & Beaumont 2002, 2003; Huismans*et al.*2005).The thermal state of the model, which also contributes to the asymmetry by determining the relative distribution of frictional–plastic and ductile regions as the model evolves (Huismans & Beaumont 2002, 2003; Lavier & Buck 2002). When extension is slow, the thermal Peclet number is small and the lithosphere cools as it extends, thereby continuously renewing the frictional–plastic layers as the model evolves. This promotes asymmetric behaviour by comparison with models that have a large thermal Peclet number (Fig. 6a).

## Analysis of simple models using the principle of minimum energy dissipation

### Approximate analytical theory

The first series of models displayed a range of behaviours that we interpreted in terms of factors that reinforced or opposed asymmetric rifting. Although this approach does provide a certain level of insight, there is no governing principle that guides the overall interpretation. In this subsection we simplify the models to consider only the most basic behaviour of a two-layer system, and take a different approach to the interpretation by demonstrating the use of the minimum rate of energy dissipation as a guide to model behaviour (Bird & Yuen 1979; Sleep *et al.* 1979; Masek & Duncan 1998). Huismans *et al.* (2005) used this approach, subject to the restrictions noted by Bird & Yuen (1979), to understand the factors that control localized asymmetric modes of deformation v. distributed symmetric modes by comparing the minimum dissipation analysis with equivalent finite-element models, which are themselves based on minimum dissipation.

We review the basic case of a two-dimensional, plane-strain system comprising a uniform, strain-softening, frictional-plastic layer overlying and bonded to a uniform linear viscous layer. The analysis predicts the rate of internal energy dissipation (henceforth termed the dissipation for brevity) during the early stages of extension for three possible modes, the pure shear (PS), symmetric plug (SP) and asymmetric plug (AP) modes (Fig. 7). Details of the analysis are given by Huismans *et al.* (2005). For this simple system there are two main contributions to the rate of internal energy dissipation. The first is from extension of the plastic layer and is the same for each of the modes except when the dissipation is reduced by strain softening of the localized shears in the SP and AP modes. The second is from the viscous layer and its size ranges from the contribution owing to pure shear, in the PS mode, to a combination of pure shear and shear of a boundary layer at the top of the viscous layer when there is differential extension between the two layers in the SP and AP modes (Fig. 7). If we consider PS as a reference mode (Fig. 7), we can then ask under what circumstances a different, lower dissipation mode exists and will be selected. It can readily be understood that any overall reduction in dissipation depends on the trade-off between the ‘gain’ of reduced plastic dissipation owing to strain softening in the SP and AP modes v. the ‘penalty’ of increased dissipation owing to development of the viscous boundary layer in the SP and AP modes (Fig. 7). By considering these trade-offs, mode selection can be predicted as a function of the layer properties and the rate of extension. The latter is important because the rate of dissipation in the plastic layer depends on the strain rate, whereas that of the viscous layer depends on the square of the strain rate, so mode selection will depend on the extension velocity, *V* (Fig. 7).

The total internal dissipation for each mode, *W*_{T}, is calculated approximately for a given choice of the model parameter values such as layer thicknesses, layer properties, the amount of strain softening and the extension velocity, and the results are plotted as a function of, for example, the viscosity, η, of the lower layer (Fig. 8). This figure shows the dissipation for the SP and AP modes. A mode transition occurs between the AP and SP modes where the curves cross at the transition viscosity, η_{T}. In this case, the SP mode becomes the lower dissipation mode for η>η_{T}. Mode transitions only occur when the dissipation curves cross; in other cases one mode systematically has a lower dissipation than another for all values of the viscosity. The total dissipation, *W*_{T}, of each mode is equal at a transition, therefore
4
where subscripts 1 and 2 successively represent the AP, SP and PS modes in our example. Solving the set of equations in equation 5 yields analytical expressions for the transition viscosities, η_{T1} and η_{T2}, which then shows how the transitions depend on the model properties (e.g. Fig. 8, equation 1, for the AP to SP transition: see Huismans *et al.* 2005, for derivation of the equations). These expressions are relatively simple if the analysis is confined to the plastic and viscous contributions, as in equation 5, but become more complex when additional contributions, such as dissipation against gravity and layer bending, are included. It can also be seen from equation 5 that these additional terms can be ignored if their contribution to both modes at a transition is approximately equal. Under this circumstance they contribute equally to both sides of the equation and can be subtracted.

The results from the approximate analysis (equation 5) demonstrate that modes are selected in the order AP, SP and PS as the viscosity of the lower layer increases and all other parameters are constant (Huismans *et al.* 2005). This order appears to be robust for the early phases of extension for parameter values appropriate for lithospheric-scale extension, and also applies to models in which the upper layer strain softens by loss of cohesion (Huismans *et al.* 2005). The effect of ignoring additional contributions to the dissipation means that the transition viscosity estimates, η_{T}, are approximate, but the errors do not affect the selection order in this particular case (Huismans *et al.* 2005). The overall results of this analysis are summarized in Figure 9, which shows how the modes are distributed in relation to the plastic and viscous controls. NSS, PSS and SS are, respectively, the non-strain-softened, partially strain-softened and strain-softened strengths of the plastic layer (the corresponding values of ϕ(ε) appear in equation 1, Fig. 8), and η_{T1} and η_{T2} are transition viscosities between AP and SP, and SP and PS modes, respectively. Arrows show how models typically evolve during strain softening of the plastic layer, when the other properties of the system, such as extension velocity, remain constant.

The mode selection and mode transitions can be understood from a physical perspective in terms of the competition between the ‘plastic gain’ and ‘viscous penalty’ contributions, defined above, to the dissipation. For example, for η>η_{T2} the mode selected is always PS (Fig. 9) because the viscous penalty incurred by the boundary layer dissipation in the SP or AP modes (Fig. 7) is always greater than the plastic gain achieved by strain softening of the localized shear zones. This result explains the perplexing observation that for this viscosity range the model will extend by pure shear even when the upper layer could develop weak localized shears. This result demonstrates that the properties of the bonded laminate are clearly different from those of its constituent layers. For η_{T1}<η<η_{T2} the strain-softened symmetric plug mode, SP2, is selected during plastic strain softening, and for η<η_{T1} the AP mode can be selected. However, as is explained below, the mode selection also depends on the way the modes evolve, that is, it is path dependent.

## Simple two-layer numerical models

The next step is to present finite-element model experiments designed to test the approximate analytical theory described in the previous section. The numerical model is the same as that used in the section on ‘Thermo-mechanical models of lithosphere extension’, but it now has the simple two-layer geometry (Fig. 10) and rheological structure (Table 2). A strain-softening frictional–plastic layer (ϕ (ε)=2°–15°) overlies and is bonded to the uniform constant linear viscous layer. The weak seed is placed at the base of the plastic layer. The results are shown after 40 km total extension at *V*=1 cm year^{−1} for viscosities, η, ranging from 10^{21} to 10^{23} Pa s (Fig. 11). At low viscosity (Fig. 11a) the AP mode is selected. Strain softening has focused the plastic deformation onto one shear and the extension is highly asymmetric. When η=10^{22} Pa s (Fig. 11b) asymmetry is suppressed and extension in the plastic layer occurs along two symmetric shears in the SP mode. At high viscosity (Fig. 11c) local deformation in the plastic layer is completely suppressed and extension of both layers occurs by the PS mode. Huismans *et al.* (2005) present a greater range of numerical experiments, but the three results (Fig. 11) serve to illustrate the basic behaviour.

### Comparison of the analytical dissipation results with the numerical models

The results of the numerical models are in qualitative agreement with the predictions of the dissipation analysis in that mode selection changes from the AP mode through SP to PS as the viscosity of the lower layer is increased (Fig. 11). The results also indicate that the transition viscosities, η_{T1} and η_{T2}, fall between 10^{21} and 10^{22} Pa s, and between 10^{22} and 10^{23} Pa s, respectively, for the parameter values (Table 2) used in the models, which is in agreement with the equivalent analytical results. The finite-element model experiments not only confirm the approximate analysis for the initial mode selection, they also demonstrate that these modes persist as the model evolves. The results also confirm our assertion that the additional contributions to the dissipation, which were ignored in the analysis but are included in the full finite-element solution, do not change the order of the mode selection, thereby providing reassurance that the approximate treatment is acceptable. The results also indicate that other modes, not considered in the analysis, are not strongly preferred for this range of viscosities. We do, however, know that modes comprising a series of plugs, not just a single plug, have levels of dissipation intermediate between the SP and PS modes. In addition, it can also be shown that the mode in which the plastic layer totally detaches from the viscous layer is not selected for the range of parameters we consider because the plastic dissipation ‘penalty’ incurred by shear on the plastic detachment is much larger than any of the possible dissipation ‘gains’ in the viscous or plastic layers for this particular mode.

The overall results can be interpreted in several ways, each of which offers a different insight into the mode selection. The first interpretation describes the mode selection in relation to increasing viscosity of the viscous layer. When the viscous dissipation is negligible, the AP mode is chosen because this provides the fastest route to minimize the plastic dissipation by strain softening only one, not two, shear zones. When viscous dissipation is larger, the mode selected will change from AP to SP at the transition viscosity η_{T2}, and to PS at η_{T2}, as explained in the previous section.

In the second interpretation the results are described in terms of the trade-off between two differential forces shown, for example, by the second equation in Figure 8. The left-hand side of the equation has dimensions of force per unit length of the model and is termed the differential ‘viscous penalty force’ incurred by deforming in the localized AP mode as opposed to the more distributed SP Mode. The right-hand side of the equation is termed the differential ‘plastic gain force’, and is equal to the reduction in the force in the plastic layer when deformed in the fully strain-softened AP mode by comparison with the PS mode. Mode transitions occur when these two forces are equal. The physical insight provided by this requirement is that the system behaviour is determined not by the intrinsic properties of either the plastic or viscous layer alone, but instead it is determined by the trade-off between the differential penalty/gain forces.

This interpretation can also expanded to consider the effect of other factors, for example, extension velocity as noted in the results of the section ‘Thermo-mechanical models of lithosphere extension’ (Fig. 4) (and also those of Huismans & Beaumont 2003) in which asymmetric, AP mode extension occurs when the rifting velocity is 0.3 cm year^{−1}, or less, but exactly the same system exhibits SP mode extension when the rifting velocity is 10 cm year^{−1}. As shown by equation 2 of Figure 8, the higher velocity modifies the viscous dissipation such that the differential viscous penalty force will be equal to the differential plastic gain force at a transition viscosity that is a factor of 33 (10/0.3 from above) smaller than that for the lower velocity, a result also implied by equation 1 of Figure 8. This modified transition viscosity is much smaller than that for the lower extension velocity. Consequently, the rapid extension model is located in the SP part of mode space (Fig. 9), whereas the model with the lower extension velocity is located in the AP part. This result demonstrates that the transition viscosities change with model properties. At even higher extension velocities the mode would change to PS. It follows from Figure 8, equation 1 that variations in any of the parameters on the right-hand side of the equation, particularly the velocity, *V*, the thickness of the plastic layer, *h*_{b}, and the initial and final strength of the plastic layer, ρ*gh*_{b}^{2} sin(ϕ_{SP})/2 and ρ*gh*_{b}^{2} sin(ϕ_{AP})/2 will modify the transition viscosity between the SP and AP modes. This result has the potential to provide insight into natural systems.

The third interpretation uses the concept of ‘dominant’ and ‘subordinate’ rheologies introduced by Huismans & Beaumont (2003). This concept is qualitative but can be made more precise by considering the differential plastic gain and viscous penalty forces. Except at a mode transition (e.g. Fig. 8, equation 2) one differential force is larger and can be interpreted to imply that the corresponding rheology ‘dominates’, whereas the other is ‘subordinate’. For example, if the viscous differential force is larger, the viscous rheology dominates and forces the model to select the mode that minimizes the viscous dissipation, and the converse. The concept provides a useful physical interpretation of the model results. The preferred mode for the plastic layer acting alone would be AP, because this is the fastest way to minimize the dissipation of this layer. This mode is not selected when the viscous rheology is dominant because the ‘dominant’ rheology forces the mode selection. The same concept applies to the velocity sensitivity of the model results described in the section on ‘Thermo-mechanical models of lithosphere extension’ (Fig. 4). At *V*=0.3 cm year^{−1}, or less, the plastic rheology dominates and the mode is therefore AP, but at *V*=10 cm year^{−1}, the viscous rheology dominates and the mode is SP.

## Systems with statistical heterogeneity

### Statistical model of inheritance from earlier tectonic deformation

We now investigate systems that are characterized by the second type of heterogeneity, statistical noise in the initial strain field. The previous two sets of models used a single weak seed to localize the deformation. As noted in the subsection on ‘Role of inheritance and heterogeneity’ earlier, this approach fails to consider that the weakest point may be ignored in favour of a distribution of inherited weak heterogeneities that link together to form even weaker shear zones. Models with statistical heterogeneities in the initial strain field are designed to represent inheritance of deformation from previous tectonic phases. They then exploit the positive feedback link between the inherited strain distribution and strain softening.

Specifically, the second invariant of the deviatoric strain, (*I*′_{2})^{1/2}, is initialized with white noise that has a truncated Gaussian distribution with a mean value ((*I*′_{2})_{mid}^{1/2}) set to half the strain at which strain softening starts ((*I*′_{2})_{min}^{1/2}), and with a maximum value just below (*I*_{2}′)^{1/2}*min*. Therefore, there is no initial strain softening in the model, but the noise preconditions parts of the model to be on the threshold of strain softening. We have also explored other types of noise, in particular fractal noise with various fractal exponents, and have found that the first order of the model behaviour is similar for a range of noise types.

### Simple two-layer models with statistical heterogeneity

The results of models that are equivalent to those of the subsection on ‘Simple two-layer numerical models’, but instead have statistical noise in the initial strain field, are shown for viscosities of the lower layer ranging from 10^{21} to 10^{23} Pa s (Fig. 12). The models all extend at *V*=1 cm year^{−1} and have exactly the same noise distribution. At low viscosity (Fig. 12a) the AP mode is selected. At intermediate viscosity (Fig. 12b) deformation is more distributed, strain softening having developed on several competing shear zones. At high viscosity (Fig. 12c) deformation in the plastic layer is distributed and characterized by a large number of shear zones. These behaviours can be interpreted using the minimum dissipation principle. The AP mode with highly asymmetric deformation can be adopted when the viscous penalty contribution to the dissipation is small (Fig. 12a), and in this case deformation is accommodated on only one plastic shear zone, thereby maximizing the plastic gain through rapid strain softening. At intermediate viscosity the mode is similar to SP (Fig. 11), but several shear zones have developed (Fig. 12b). Evidently, the total dissipation is best minimized by the activation of several plastic shears even though some will strain soften faster than others. In this case the major advantage is the ‘gain’ of reduced viscous dissipation by comparison with the AP and SP modes. Similar reasoning suggests the SP mode with a single plug should be selected for viscosities somewhere in the range 10^{21} Pa s < η<10^{22} Pa s. At high viscosity, 10^{23} Pa s, the model (Fig. 12c) is dominated by the viscous rheology, at least in a relative sense. Some plastic gain has been achieved by strain softening on the plastic shear zones, but it is at the expense that many shears had to be created. The major gain is in the corresponding reduction in the viscous penalty because there are now only short, low strain rate, segments of the boundary layer flow between the multiple plastic shear zones, not long segments that accommodate major offsets as in the AP and SP modes.

The overall results indicate that this heterogeneous noise, which is below the strain-softening threshold, does not dictate the behaviour of the system and that deformation similar to the AP, SP and PS modes can occur depending on the viscosity of the lower layer. The difference is that modes with multiple plastic shear zones are now easily excited for the viscosity range η_{T1}<η<η_{T2} because the inherited strain heterogeneity places many regions of the model close to the threshold for strain softening. These new modes modify Figure 9 by occupying part of the mode space previously taken by the SP and PS modes, such that a new transition viscosity, to the multiple plastic shear modes, occurs between η_{T1} and η_{T2}, and the value of η_{T2} is increased. The multiple plastic shear modes will be characterized by an increasing number of shears with increasing viscosity of the underlying layer. This behaviour occurs because the viscous penalty is reduced by reducing the spacing between the plastic shears, as noted above, and in the limit of high viscosity, the new η_{T2}, the overall deformation will approach the PS mode.

### Lithosphere-scale extension with statistical crustal heterogeneity

In this section we investigate the role of statistical heterogeneity in the initial strain field in lithosphere-scale models equivalent to those of the section on ‘Thermo-mechanical models of lithosphere extension’ to determine whether mode selection will be modified. We assume heterogeneity, due to earlier phases of deformation, is thermally annealed in the mantle where temperatures are high and, therefore, only apply statistical heterogeneity to the crust. Note that strain softening does not occur in the viscous (ductile) regime so that models with a strong (plastic) lower crust will experience strain softening in this region but those with a weak, ductile lower crust will not. For this series of models (*I*′_{2})^{1/2} is initialized with white noise that has a Gaussian distribution with a mean value ((*I*′_{2})_{mid}^{1/2}) set to the strain at which strain softening starts ((*I*′_{2})_{min}^{1/2}), and with the same Gaussian distribution. The other properties of these models are the same as the equivalent ones in the section on ‘Thermo-mechanical models of lithosphere extension’ except that the extension velocity is 0.5 cm year^{−1}.

When the lower crust is strong (Fig. 13) extension is at first pure shear. When the incremental strain augments the initial statistical strain to the threshold for strain softening, strain softening occurs in the frictional-plastic crust and the SP mode is selected. The position of the plug depends on the particular realization of the white noise, a behaviour that is equivalent to randomly placing a weak seed. Localized deformation in the crust causes diffuse projections of the PS mode shears to develop in the upper mantle lithosphere and these forced shear zones eventually localize by strain softening in the upper frictional–plastic mantle. Feedback of strain softening results in the development of an asymmetric shear zone that penetrates most of the lithosphere (Fig. 13a) which facilitates exhumation of the upper mantle lithosphere (Fig. 13b) via the AP mode. During later stages, when thermal advection weakens the lithosphere, viscous necking becomes the dominant style of deformation (Fig. 13c). This model behaviour is very similar to its single weak seed equivalent (Figs 4b & 5b). The difference in the extension velocity, 0.5 v. 0.3 cm year^{−1}, has only a minor effect on results. A more important implication is that models with this particular set of properties behave in a similar manner and are not particularly sensitive to the difference between inherited statistical heterogeneity and a deterministic weak seed. The result is also consistent with the simple two-layer model with statistical heterogeneity and a low viscosity lower layer (Fig. 12a). Evidently, in all of these examples, which are essentially two-layer cases, the plastic gain achieved by selecting the AP mode outweighs the viscous penalty for low values of the effective mantle viscosity. For higher effective viscosities of the mantle lithosphere we may expect more distributed modes of deformation like those seen in Fig. 12b, c.

The style of deformation for the equivalent model with the viscously weak lower crust (Fig. 14) is markedly different, both from the model with strong lower crust and, more significantly, from the equivalent model with a single weak seed (compare Figs 5a and 14). At first deformation occurs in an overall ‘PS-type’ mode (Fig. 14a). Sets of parallel and conjugate frictional shears form in the frictional upper crust, but strain localization does not occur in the upper mantle lithosphere which has no initial statistical heterogeneity. After approximately 160 km of extension, localization in the crustal layer feeds back into the frictional–plastic upper mantle lithosphere and seeds a localized shear zone (Fig. 14b). The remaining evolution is similar to the equivalent model with a strong lower crust except that the final necking is more symmetrical and there is somewhat less exhumation of the mantle (Fig. 14c). The significance of this model is that it displays a natural progression from an early distributed wide rift mode to a late localized narrow rift mode. Although the equivalent model with the single weak seed also exhibits decoupling and differential shear between the crust and uppermost mantle (Fig. 5a), the width of the zone of distributed crustal deformation is much narrower and tends to increase as rifting progresses.

## Implications for the Iberia–Newfoundland conjugate margin system

### First-order characteristics of the Iberia–Newfoundland conjugate margin system

In this subsection we explore whether what has been learned from the models above can be applied to extension of the Iberia–Newfoundland conjugate margin system. One of the main characteristics of the Iberia–Newfoundland conjugate margin is the high degree of variability in the geometry of the margins and the location of extension with position along the margins (Figs 15 & 16). This variability may result from the competition between local inherited structures and underlying physical processes that control the mode of extension. We first attempt to identify the general features of the system because these may lead to insight into the fundamental underlying processes responsible for its formation. Secondly, we address the potential causes of the variability.

Cross-sections (Keen *et al.* 1987*a*, *b*; Torne *et al.* 1994; Dean *et al.* 2000; Funck *et al.* 2003; Pérez-Gussinyé *et al.* 2003) shown restored to approximate the configuration at the time of onset of sea-floor spreading (Fig. 16) illustrate the strike variability over a distance of only 100 km. We regard the following to be key features of the Iberia–Newfoundland conjugate margin system. (1) Extension and rifting occurred during a prolonged period of approximately 80 Ma with at least three phases of rift activity (Mauffret & Montadert 1988; Tankard & Balkwill 1989; Tankard *et al.* 1989; Murillas *et al.* 1990; Driscoll *et al.* 1995; Stapel *et al.* 1996). (2) Initial rifting occurred in a wide rift mode over an extended area that was at least 400 km wide (Keen *et al.* 1987*a*, *b*; Tankard & Welsink 1989; Manatschal 2004). (3) The system became a narrow rift late in its evolution and continental mantle lithosphere was exhumed during the terminal phase (Boillot *et al.* 1988; Louden & Chian *et al.* 1999; Dean *et al.* 2000). (4) The system is a non-volcanic rift zone, the late stages of rifting and early sea-floor spreading being remarkably a-magmatic (Louden & Chian 1999). (5) Extension occurred at very low average velocities, an important point that we substantiate below. (6) Localized extension during the first two phases of rifting was limited to the crust with only minor thinning of the mantle lithosphere as evidenced by limited post-rift subsidence of the early rift basins (Stapel *et al.* 1996). (7) Asymmetries in crustal thickness, distribution of exhumed mantle lithosphere (Whitmarsh *et al.* 2001) and spatial variability of mineral chemistry of the exhumed mantle lithosphere (Muentener 2004) suggest that late-stage rifting was accommodated by an asymmetric shear that cut down into the mantle lithosphere (see also Whitmarsh *et al.* 2001; Manatschal 2004).

Timing, duration and associated extension of the rift phases are not well constrained. We have constructed bounds on extension velocity and duration of each of the rift phases using various permissible rift-phase reconstructions and estimates of the total amount of crustal extension. The total crustal extension is estimated using the northern conjugate margin transect (Funck *et al.* 2003; Pérez-Gussinyé *et al.* 2003) as a reference. Total extension of this part of the margin, prior to the onset of sea-floor spreading, ranges from 180 to 220 km, based on the cross-sectional area balance between an estimated original crustal thickness of 32 km and the present-day thickness distribution of the continental crust. This range overlaps with that based on a similar estimate of 100–200 km derived from cross-sections in the central part of the Iberia–Newfoundland conjugate margin system (Minshull *et al.* 2001), although it is clear their estimate has a large uncertainty. This total extension must be apportioned among the different phases.

The timing of the rift phases can be determined approximately based on subsidence, from well data and seismic profiles (Tankard & Balkwill 1989), as a proxy for extension:

- Phase I:
- 218–205 Ma – minor relief formation in external basins;
- Phase II:
- 165–142 Ma – major relief formation in external basins (e.g. Jeanne d'Arc and associated basins, Galicia Interior Basin, Lusitanian Basin);
- Phase III:
- 142–122 Ma – Narrow rifting and breakup of the crust;
- Phase IV:
- 122–116 Ma – Asymmetric exhumation of the mantle lithosphere.

Previous estimates of the extension velocity are of the order of 1.5 cm year^{−1} full spreading velocity (Minshull *et al.* 2001). The approximate timing of rift phases and the value of 200 km total extension given above can be used to make simple estimates that give upper and lower bounds for the average extension velocity of the system. If the extension velocity was the same throughout all of the phases the rate was approximately 0.3 cm year^{−1} full spreading rate. If all of the extension occurred during the latter two phases the estimated velocity would be approximately 0.75 cm year^{−1}. We regard this estimate to be an unreasonably large upper bound because phases I and II achieved significant normal faulting and basin subsidence. We therefore regard 0.5 cm year^{−1} to be a representative upper-bound estimate of the average extension velocity during the phases when extension occurred.

In summary, from the discussion above, a successful model needs to meet the following constraints:

an initially wide rifting mode with extension/subsidence focused in a number of distinct, widely distributed, fault-bounded basins, with structure that varies with position along the system;

Initial extension in these basins was mostly limited to the crust;

Abandonment of the external basins in favour of late-stage narrow, possibly asymmetric, rifting;

exhumation of the mantle lithosphere during the final phases of breakup;

An overall slow extension rate probably close to 0.5 cm year

^{−1}.

### Towards a first-order model with the characteristics of the Iberia–Newfoundland conjugate margin system

In the context of the models described in this paper, the progression from a wide to narrow rift can be met by models that inherit heterogeneous crustal noise and in which the crustal noise develops into a wide region of extension characterized by weak upper crustal shears and faults that bound distinct basins. Horizontal decoupling and shear in the lower or mid crust appear to be necessary, in addition to the inherited heterogeneity, (cf. Fig. 14) because the wide rifting mode does not develop when the lower crust is strong and precludes decoupling (cf. Fig. 13). The progression from a wide to a narrow rift occurs naturally in this type of model if the uppermost mantle lithosphere strainsoftens, leading to localization of the deformation (cf. Figs 13 & 14). It also appears that models with a strong, coupled, lower crust show the greatest potential to exhume mantle lithosphere to the surface during the narrow rifting mode (e.g. Figs 5b & 13c) and that this exhumation is commonly achieved by asymmetric extension (cf. Fig. 4b, c). These behaviours are also compatible with slow rifting velocities that favour AP and SP modes over distributed pure shear, as described in the section on ‘Analysis of simple models using the principle of minimum energy dissipation’.

The model shown in Figure 14 has many of the required general characteristics displayed by the Newfoundland–Iberia conjugate margin system, as listed above. However, it is a model with a weak lower crust and, therefore, fails to exhume the mantle to the surface during the late narrow rifting phase. To date, we have been unable to create a self-consistent model that produces both initial wide rifting and later stage cooling that strengthens the viscous lower crust sufficiently to create the necessary conditions for very narrow AP or SP rifting and mantle exhumation. We suspect this inability occurs because we consider models in which extension is uniform and continuous, and, therefore, there are no tectonically quiescent periods during which significant cooling can occur.

We demonstrate the effect of cooling and strengthening of the extended crust using a composite model. The rift evolution starts in a model with inherited strain heterogeneity and a weak middle and lower crust. We then simulate the effects of cooling and late-stage embrittlement of the crust by suppressing the viscous flow in the crust during Phase 2. This composite behaviour is obtained by increasing crustal viscosity by a factor 100 during Phase 2 to achieve a crust that is totally in the frictional–plastic regime. The extension velocity is 0.5 cm year^{−1}. During the 24 Ma-long Phase 1, 120 km of extension is accommodated by localization on multiple sets of frictional shear zones in the upper crust (Fig. 17a). Deformation of the mantle lithosphere occurs in a distributed mode. During Phase 2, localization in the crustal layer feeds back into the upper mantle lithosphere where it causes localization in the frictional–plastic uppermost mantle. Deformation is asymmetric (Fig. 17b) and mantle lithosphere is exhumed during the final phases of breakup (Fig. 17c).

To first order, this model result is compatible with the constraints compiled for the Newfoundland–Iberia natural system (Fig. 16), particularly when the significant variability among cross sections located close together is taken into account (Fig. 16a–c). It also provides a potential explanation for this variability in the geometry of the margins, in particular the initial irregular distribution and later abandonment of early rift basins that flank the margins if these are the consequence of deformation related to inherited crustal heterogeneity. We anticipate that an equivalent fully three-dimensional (3D) model, with inherited initial heterogeneous strain, would develop sets of anastamosing rift basins during the early phase of rifting and that the structure of these basins would vary significantly along strike of the system. These preliminary results lead us to the proposition that the early stages of rifting in margins of this type are disorganized, reflecting inherited lithospheric characteristics. Only later in the rifting does the system overcome this inheritance and become organized in the form of a narrow rift.

## Summary and conclusions

This article has focused on models of extensional processes with application to lithospheric extension during rifting and the formation of rifted continental margins. Even though rifts and rifted margins have been studied intensively for more than 30 years, a fundamental understanding of the factors that control the extensional process, the distribution of extension and the associated rift and rift margin geometry is still lacking. The purpose of the model studies reported here and in associated articles (Huismans & Beaumont 2002, 2003; Huismans *et al.* 2005) is to test a limited range of conceptual models that combine distributed viscous–plastic flows with feedback mechanisms, specifically plastic strain softening, that lead to localized deformation. This localization allows the models to represent faults and shear zones approximately. In our earlier research, the development of localized shears was in part deterministic and was triggered by *a priori* inclusion of small weak regions, ‘weak seeds’, in the models. In this paper we have reviewed those results and expanded the research to include the effect of distributions of weak heterogeneities that take the form of white noise in the initial strain field and can be interpreted to represent strain, or in a wider sense damage, that is inherited. The effects of the viscous‐plastic layered rheology combined with the inherited initial strain lead to more general representations of the type of deformation that can occur during lithospheric extension. We believe this approach, which combines a continuum model based on plastic and viscous rheologies with strain softening that results from a statistical distribution of inherited heterogeneity, provides a potential explanation for the general features of the Newfoundland–Iberia conjugate margins. These margins appear to have had an early phase of distributed disorganized crustal extension, and only later became organized into a narrow rift that exhumed continental mantle lithosphere during terminal extension prior to the onset of sea-floor spreading.

## Acknowledgments

We wish to thank the organizers of the IMEDL workshop for the invitation to participate. This research was funded through an ACOA–Atlantic Innovation Fund contract, and an IBM-Shared University Research grant. C. Beaumont was funded by the Canada Research Chair in Geodynamics. C. Beaumont also acknowledges support through the Inco Fellowship of the Canadian Institute for Advanced Research. Numerical calculations used software developed by P. Fullsack. We thank R. Buck, L. Lavier and an anonymous reviewer for the constructive and thorough reviews.

- © The Geological Society of London 2007