Geodynamic Models

32 minutes read

McKenzie Stretching Model

The McKenzie model is the simplest stretching and cooling model, the easiest to understand, and the basis for more sophisticated models. It is useful to discuss the assumptions, limitations, and advantages of the McKenzie model first, before moving on to more complex models. We focus on the physical approximations and predictions of the McKenzie model and their implications for basin evolution and analysis. For a more complete discussion of the model derivations, refer to McKenzie (1978). Table 1 defines the relevant parameters.

Parameter Value Used in Examples

Density of Water	ρ_w	1.00 gm/cm³
Density of Crust at 0º C	ρ_Co	2.80 gm/cm³
Density of Crust	ρ_c	– calculated – gm/cm³
Density of Mantle Lithosphere	ρ₁	– calculated – gm/cm³
Density of Mantle Lithosphere	ρ_1o	3.33 gm/cm³ at 0º C
Density of Mantle Lithosphere	ρ_m	3.18 gm/cm³ at the Solidus
Thickness of Lithosphere a 125 km
Initial Thickness of Crust	t_c	-variable – (km)
Lithospheric Stretching Factor	β	– variable – (dimensionless)
Surface Temperature of	T_o	0º C
Lithosphere
Solidus Temperature of	T_m	1333º C
Lithosphere
Volume Coefficient of Thermal	α	3.28⋅10⁻⁵ ^∘C⁻¹
Expansion
Thermal Time Constant	τ	62.8 My
Thermal Diffusivity	κ	8⋅10−7 m²/s
Equilibrium Surface Heat Flux	κTm/a	0.8 μcalcm⁻²s⁻¹

Table 1: Parameters and their values used in the discussion of cooling and stretching models.

The original McKenzie model (Figure 1)

The original McKenzie model, Geodynamic Models — **FIGURE 1** The original McKenzie model

Figure 2,

after rifting, the initial temperature distribution is altered only by a raised lithosphere-asthenosphere boundary, Geodynamic Models — **FIGURE 2** After rifting, the initial temperature distribution is altered only by a raised lithosphere/asthenosphere boundary

and Figure 3) assumes that the ductile and brittle portions of the lithosphere stretch equally.

the lithosphere gradually cools and returns to its initial thickness and geothermal gradient, Geodynamic Models — **FIGURE 3** The lithosphere gradually cools and returns to its initial thickness and geothermal gradient

The McKenzie model assumes no rheological differentiation or inhomogeneity within the lithosphere. The mantle lithosphere and crust stretch by a factor β and, due to conservation of volume, thin by a factor 1β with asthenospheric material rising to fill the created space. The McKenzie model ignores rifting dynamics by assuming instantaneous stretching and thinning, followed by post-rift cooling and contraction with one-dimensional conductive heat flow. Immediately after rifting, the initial temperature distribution is altered only by a raised lithosphere/asthenosphere boundary (Figure 2). With the boundary closer to the earth’s surface, the geothermal gradient in the lithosphere increases. Subsequently, the lithosphere gradually cools and returns to its initial thickness and geothermal gradient. (Figure 3). The McKenzie model estimates surface elevation changes caused by Airy isostatic compensation of thermal loads. These thermal loads are changes in the vertical density structure induced by temperature changes caused by stretching and then cooling of the lithosphere.

We estimate elevation changes in two stages. Initial changes are due to one-step instantaneous rifting, termed initial subsidence (Figure 2; model time equals zero). Subsequent changes through time are due to conductive cooling, termed thermal subsidence (Figure 3; model time greater than zero). The total subsidence, S(t), termed tectonic subsidence , is the sum of the initial subsidence, S_i, and the thermal subsidence, S_T(t):

$S(t) = S_i + S_T(t)$ $\quad$ (1)

For each model stage, we calculate densities from the vertical temperature structure. The density of crust and mantle lithosphere with depth and time, ρ_c(z,t) and ρl(z,t), are related to temperature, T(z,t), by α, a coefficient of thermal expansion:

$\rho_c(z,\,t)=\rho_{c_{0}} \left[ 1 - \alpha (T(z,\, t) - T_0) \right]$ $\quad$ (2)

$\rho_1(z,\,t)=\rho_{1_{0}} \left[ 1 - \alpha (T(z,\, t) - T_0) \right]$ $\quad$ (3)

We then calculate surface elevation changes through time by applying Airy compensation to the changing density structure, ρ_c(z,t) and ρ₁(z,t). In the following derivations, we assume T₀ is 0 degrees C and consequently drops out of Equations 2 and 3.

Using this approach, initial subsidence, Si, is given by

$S_i = \dfrac{a\, \left[ (\rho_m - \rho_{c_{0}})\dfrac{t_c}{a}\left( 1 - \dfrac{\alpha\, T_m\, t_c}{2a} \right) - \dfrac{\rho_m\, \alpha \, T_m}{2} \right] }{\rho_m \left( 1- \alpha \, T_m \right) - \rho _w} \left( 1 - \dfrac{1}{\beta} \right)$ $\quad$ (4)

In this derivation, water is the replacement mass , the mass that fills the topographic depression due to subsidence. Si is positive for negative surface elevation changes (subsidence) and negative for positive surface elevation changes (uplift).

One way to understand how subsidence or uplift occurs as a response to lithospheric stretching is to consider the changes in relative buoyancy that occur when lithosphere thins and the resulting lithospheric volume is replaced by asthenosphere. As shown in Figure 4,

relative-buoyancy that occur when lithosphere thins and the resulting lithospheric volume is replaced by asthenosphere, Geodynamic Models — **FIGURE 4** Relative-buoyancy that occur when lithosphere thins and the resulting lithospheric volume is replaced by asthenosphere

we can rewrite Equation 4, incorporating the hydrostatic balance of an incremental load due to an extended lithosphere supported by a fluid asthenosphere (modified for water loading from Beaumont et al., 1982):

$S_i = \left[ \dfrac{\rho _m\, a - \rho_{cav}\, t_c - \rho_{lav}(a - t_c)}{\rho_m - \rho_w} \right] \left( 1 - \dfrac{1}{\beta} \right)$ $\quad$ (5)

where

$\rho_{m} = \rho_0(1-\alpha\, T_m)$ , the density of the lithosphere at the solidus

$\rho_{cav} = \rho_{c_{0}}(1-\alpha \, T_{cav})$ , the average density of the crust

$\rho_{lav} = \rho_0(1-\alpha \, T_{lav})$ , the average density of the mantle lithosphere

$T_{cav}= \tfrac{T_m\, t_c}{2a}$ , the average temperature of the crust; and

$T_{lav} = \tfrac{T_m}{2}(1+\tfrac{t_c}{a})$ , the average temperature of the mantle lithosphere

The parameters with the highest variability in the McKenzie model are the initial crustal and lithospheric thicknesses and the degree of stretching. For a specific crustal and lithospheric configuration, the McKenzie model predicts initial subsidence or uplift primarily as a function of the degree of stretching, β. Initial subsidence increases with increased degree of stretching, asymptotically approaching a maximum value defined by the bracketed expression of Equation 5 (Figure 5 and Figure 6).

The parameters with the highest variability in the McKenzie model, Geodynamic Models — **FIGURE 5** The parameters with the highest variability in the McKenzie model

FIGURE 5

As β increases, more light crust is replaced by heavier asthenosphere, thereby inducing more subsidence.

As initial crustal thickness, t_c, decreases, the value of the bracketed expression in Equation 4 decreases, eventually becoming negative. In this situation, because less light crust is replaced by asthenosphere, either initial subsidence is reduced or there is initial uplift. Therefore, Equation 5 predicts initial uplift for a thinner crust, and subsidence for a thicker crust. For typical parameter values, uplift will occur for initial crustal thicknesses of less than approximately 17 km.

After stretching and thinning increase the temperature gradient in the lithosphere, the lithosphere cools conductively to the initial geothermal gradient (Figure 3 and Figure 8, Theoretical heat flux at the surface of the earth as a function of β from the McKenzie model [equation 7].

Theoretical heat flux at the surface of the earth as a function of β from the McKenzie model, Geodynamic Models — **FIGURE 8** Theoretical heat flux at the surface of the earth as a function of β from the McKenzie model

Values for β range from 1.25 to 100. Initial pulse of heat decays rapidly within the first 20 my after stretching. Initial surface heat flow is greater for larger values of β.) We estimate the temperature, T(z,t), and heat flow histories, F(z,t), during cooling by applying a one-dimensional conduction model:

$T(z,t) = T_m\left( \dfrac{z}{a} \right) + T_m \left( \dfrac{2}{\pi} \right) \displaystyle \sum_{n=1}^{\infty} \left[ \dfrac{\beta}{n^2\, \pi} \sin \left( \dfrac{n\, \pi}{\beta} \right) \right] \sin \left( \dfrac{n\, \pi\, z}{a} \right)e ^{\tfrac{-n^2\, t}{\tau}}$ $\quad$ (6)

$F(z=surface, t) = \dfrac{\kappa \, T_m}{a} \left[ 1 + 2 \displaystyle \sum_{n=1}^{\infty} \left( \dfrac{\beta}{n^2 \pi} \sin \dfrac{n\, \pi}{\beta} \right) e ^{\tfrac{-n^2\, t}{\tau}} \right]$ $\quad$ (7)

where

$\tau = \dfrac{a^2}{\pi^2\kappa}$

Because these expressions are modified exponentials, they predict rapid decay of heat flow and equilibration to the original geothermal gradient (Figure 7), Theoretical temperature gradients resulting from thermal decay to an equilibrium steady state after instantaneous stretching from the McKenzie model [equation 6].

Theoretical temperature gradients resulting from thermal decay to an equilibrium steady state after instantaneous stretching from the McKenzie model, Geodynamic Models — **FIGURE 7** Theoretical temperature gradients resulting from thermal decay to an equilibrium steady state after instantaneous stretching from the McKenzie model

The equilibrium thermal gradient equals the gradient before instantaneous stretching. Immediately after instantaneous stretching, the solidus temperature of the mantle is raised, thereby increasing the geothermal gradient. As the lithosphere cools conductively, the geothermal gradient approaches the original equilibrium gradient. Because thermal decay is initially rapid, but decreases approximately exponentially with time, the curves approach the equilibrium gradient most rapidly just after stretching.) and (Figure 8). There is an initial pulse of heat that decays rapidly early in the history of the basin. For common parameter values, the McKenzie model predicts most heat loss within the first 20 My after rifting. Larger degrees of stretching bring hot asthenospheric material nearer to the surface, thereby inducing a larger heating event. Consequently, Equations 6 and 7 predict a larger initial geothermal gradient and surface heat flow for larger values of β. For smaller degrees of stretching, heat flow early in the history of the basin is only slightly higher than later heat flow.

Similar to initial subsidence, we can estimate thermal subsidence through time, S_T(t), by calculating a vertical density structure through time from the temperature history using Equations 2, 3 and 6. We determine thermal subsidence by applying an Airy compensation model to the calculated density structure. In this derivation, the depression caused by thermal subsidence is filled with water:

$S_T(t) = \dfrac{a\, \rho_0\, \alpha \, T_m}{\rho_0 - \rho_w}\left \{ \dfrac{4}{\pi^2} \displaystyle \sum _{m=0}^{\infty }\dfrac{1}{(2\, m + 1)^2}\left [ \dfrac{\beta}{(2\, m + 1)\pi} \sin \dfrac{(2\, m +1)\pi}{\beta} \right ]\left [ 1-e^{-(2\, m + 1 )^2 \dfrac{t}{\tau}} \right ]\right \}$ $\quad$ (8)

Although Equation 8 is mathematically more complicated than Equation 5, Equation 8 expresses a physical approximation that is analogous to that expressed in Equation 5. Thermal subsidence results from the buoyant balance due to a cooling and thickening lithosphere replacing a less dense underlying fluid asthenosphere. Equation 8 predicts that the form of thermal subsidence is approximately exponential, with subsidence rates highest early in the history of the basin and gradually decreasing through time (Figure 5 and Figure 6, Upper graph shows theoretical initial and thermal subsidence as a function of degree of stretching, β. Subsidence values are referenced to the elevation when thermal subsidence initiates. To determine the initial subsidence for a value for S_i on the “Subsidence Magnitude” scale bar to the left. For example, for β equals 2, initial subsidence equals approximately 1.5 km. The lower graph shows initial and thermal subsidence [i.e., tectonic subsidence] in an absolute reference frame for β equals 2. Although the initial subsidence that occurs during rifting is considered instantaneous in the McKenzie model, as illustrated in the upper graph, the lower graph illustrates linear subsidence during a finite time span.). Because the heating event is greatest for higher degrees of stretching, thermal subsidence magnitudes are also greater for greater values of β . This form of rapid early subsidence with high heat flow that decays quickly is an important influence on the stratigraphic architecture and maturation histories within basins formed by lithospheric extension.

To apply the McKenzie model to a specific basin, we must estimate or assume the initial lithospheric thickness, time of initiation of thermal subsidence, and degree of stretching, β . One approach for estimating these parameters is through pseudo-inverse procedures (Lerche, 1990), or “curve matching,” in which we use observations or analyses of observations to find a best fit for each parameter. For example, we might use backstrip analyses to estimate subsidence. Once we obtain these values, we use the McKenzie model to estimate subsidence and heat flow for times and locations beyond the available data. Commonly, we use only subsidence analyses for estimating model values because maturation levels are much less sensitive measures of heat flow than backstrip analyses are of subsidence, and the kinetics of organic maturation parameters are poorly known. After we determine the model parameters from subsidence curve matching, we use the model to predict a heat flow history.

We can constrain estimates of β by assuming uniform stretching and that the initial crustal thicknesses in the basin were equal to the present unstretched crustal thickness on the adjacent continental margin. We can calculate β directly if we have seismic or gravity measurements of crustal thicknesses within the basin and on the unstretched margin.

The inverse method may be difficult to apply if there have been large amounts of stretching, or if available data are only from periods late in the history of the basin. We can estimate model parameters from data only if differences in curve shapes and magnitudes of model predictions are larger than error estimates. Figure 5, Figure 6, and Figure 8 illustrate that values for heat flow and subsidence converge asymptotically for large values of β, thus making it more difficult to invert for parameters as β increases. Similarly, heat flow estimates for all values of β converge with time, making the shapes of subsidence curves very similar late in the basin history. Consequently, if there are no data for estimating the shape of early subsidence, we may not be able to may to distinguish between the magnitudes of initial and thermal subsidence, making it difficult to estimate model parameters.