To model the extension of continental lithosphere, we used the finite-element method (FEM) originally employed by Dunbar and Sawyer (1988, 1989) and also by Harry and Sawyer (1992a, 1992b). (A more complete treatment of the modeling method can be found in these papers.) In this method, a vertical slice of the lithosphere was represented as a two-dimensional thermomechanical continuum, approximated by a "mesh" of elements.
The deformation of the model was governed by two different empirically determined rheological laws. The yield stress of each element was computed using Byerlee's (1978) law frictional sliding and also using a power law viscous creep of the form
The value for yield stress actually used for each element at each time step was the lower of the values predicted by these two laws, because deformation most likely occurs by the mechanism that requires the least stress. The temperature, pressure, strain rate, and lithology all influence the deformation mechanism for which the yield stress is least.
Brittle faulting in the upper crust was approximated by ideal plasticity, and no deformation by discrete faulting occurred in the model. Thus, although the formation of individual fault-controlled basins was not computed, the overall extension in the upper crust was estimated. The FEM included neither sedimentation nor any magmatic process such as partial melting. The FEM included, however, the effects of buoyancy forces, the model's own flexural strength, and the weight of any overlying seawater. The driving forces of plate tectonics, and continental lithospheric extension in particular, were approximated by constant-velocity boundary conditions on the sides of the model (Keen, 1985; Kusznir and Ziegler, 1992; Ziegler, 1992).
The FEM included heat generation in the crust; the amount of heat generated per unit volume decreased exponentially with depth. The top of the model was kept at 0°C throughout the model run. At the bottom of the model lithosphere, a constant heat flow (e.g., Braun and Beaumont, 1989a, 1989b) was maintained. (The method used to calculate the heat-flow values is described below.) We chose this boundary condition because we were interested in modeling the effect of a "resting phase" on the pattern of extension; the conductive cooling of extended lithosphere during this phase would have occurred too slowly, if at all, if a constant-temperature bottom boundary condition (e.g., Braun and Beaumont, 1987) had been chosen.
The bottom-temperature boundary condition was determined in the following manner. The unextended initial model was allowed to equilibrate thermally with a temperature of 0°C at the top and 1333°C at the bottom of the model. The vertical temperature gradient across the bottom row of elements was multiplied by the mantle conductivity to arrive at the vertical heat flow. Because the lithospheric thickness—in effect, the distance separating the 0° and 1333°C isotherms—could vary horizontally across the model, so did the thermal gradient; hence, the heat flow could vary laterally across the model as well. These heat-flow values were used as the bottom boundary condition for the remainder of the model run.
It must be noted that the constant heat-flow boundary condition will probably permit the lithosphere to cool too quickly. This boundary condition does not allow upwelling asthenosphere to become part of the lithosphere, an event that we think may have occurred during rifting of the Newfoundland-Iberian Margins. The modeling routine employed here offered only two choices of bottom boundary condition—constant temperature and constant heat flow—and modifying the code to model the thermal regime more accurately was outside the scope of the research described here. In cases of prolonged or interrupted rifting, we believe that the choice of constant heat flow is clearly the more accurate of the two options.
The modeling routine also lacks the blanketing effects of sedimentation. If sediments, which generally have low conductivity, had been included in the thermal model, the cooling of the lithosphere would have occurred at a slower rate, and thus the patterns of rifting may have been different. We are unsure of the implications of this model feature on our results; nevertheless, the conclusions of this paper still hold true for the thermal regimes used.
The continental lithosphere was represented by three rock types and their associated rheological properties (Table 1). The model mantle was wet Aheim dunite (Chopra and Paterson, 1981). Two lithologies were used to model continental crust: wet quartz diorite (Hansen and Carter, 1982) represented typical continental crust and wet granite (Hansen and Carter, 1983) represented weakened continental crust. These rheologies had been used previously by Harry and Sawyer (1992a) to model the rifting lithosphere. We did not investigate the effect of using different rheologies on multiphase rifting. The thermal properties listed in Table 2 were also used in all models.