# Continental extension#

*This section was contributed by John Naliboff, Anne Glerum, and Valentina
Magni*

In the crustal deformation examples above, the viscosity depends solely on the Drucker Prager yield criterion defined by the cohesion and internal friction angle. While this approximation works reasonably well for the uppermost crust, deeper portions of the lithosphere may undergo either brittle or viscous deformation, with the latter depending on a combination of composition, temperature, pressure and strain-rate. In effect, a combination of the Drucker-Prager and Diffusion dislocation material models is required. The visco-plastic material model is designed to take into account both brittle (plastic) and non-linear viscous deformation, thus providing a template for modeling complex lithospheric processes. Such a material model can be used in using the following set of input parameters:

```
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
end
end
```

This cookbook provides one such example where the continental lithosphere
undergoes extension. Notably, the model design follows that of numerous
previously published continental extension studies ([Brune *et al.*, 2014, Huismans and Beaumont, 2011, Naliboff and Buiter, 2015] and references therein)

## Continental Extension#

The 2D Cartesian model spans 200 (\(x\)) by 100 (\(y\)) km and has a finite element grid with 1.25 and 2.5 km grid spacing, respectively, above and below 50 km depth. This variation in grid spacing is achieved with a single initial adaptive refinement step using the minimum refinement function strategy. Unlike the crustal deformation cookbook (see Crustal deformation), the mesh is not refined with time.

```
subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 10
set Y repetitions = 5
set X extent = 200e3
set Y extent = 100e3
end
end
subsection Mesh refinement
set Initial adaptive refinement = 1
set Initial global refinement = 3
set Time steps between mesh refinement = 0
set Strategy = minimum refinement function
subsection Minimum refinement function
set Coordinate system = cartesian
set Variable names = x,y
set Function expression = if ( y>=50e3 && x>=40.e3 && x<=160.e3, 4, 3)
end
end
```

Similar to the crustal deformation examples above, this model contains a free surface. However, in this example the free surface is advected using the full velocity (e.g., normal projection) rather than only the vertical component. As this projection can lead to significant surface mesh deformation and associated solver convergence issues, diffusion is applied to the free surface at each time step. Deformation is driven by constant horizontal (\(x\)-component) velocities (0.25 cm/yr) on the side boundaries (\(y\)-velocity component unconstrained), while the bottom boundary has vertical inflow to balance the lateral outflow. The top, and bottom boundaries have fixed temperatures, while the sides are insulating. The bottom boundary is also assigned a fixed composition, while the top and sides are unconstrained.

```
subsection Boundary velocity model
set Prescribed velocity boundary indicators = left x: function, right x:function, bottom y:function
subsection Function
set Variable names = x,y
set Function constants = v=0.0025, w=200.e3, d=100.e3
set Function expression = if (x < w/2 , -v, v) ; v*2*d/w
end
end
subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface, top: diffusion
subsection Free surface
set Surface velocity projection = normal
end
subsection Diffusion
set Hillslope transport coefficient = 1.e-8
end
end
subsection Boundary composition model
set Fixed composition boundary indicators = bottom
set List of model names = initial composition
end
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 1613
set Top temperature = 273
end
end
```

Sections of the lithosphere with distinct properties are represented by
compositional fields for the upper crust (20 km thick), lower crust (20 km
thick) and mantle lithosphere (60 km thick). Material (viscous flow law
parameters, cohesion, internal friction angle) and thermodynamic properties
for each compositional field are based largely on previous numerical studies.
Dislocation creep viscous flow parameters are taken from published deformation
experiments for wet quartzite ([Rutter and Brodie, 2004]), wet anorthite ([Rybacki *et al.*, 2006]) and dry olivine ([Hirth and Kohlstedt, 2004]). Additional
compositional fields are used to track plastic strain and the non-initial
plastic strain, with the latter value tracking the same quantity as the
plastic strain absent the initial plastic strain values. As discussed further
on, the plastic strain is used to soften (e.g., reduce) the friction and
cohesion through time based on user-specified bounds and magnitudes. The
initial randomized values of plastic strain in the model center localize
distributed deformation in this region.

```
subsection Compositional fields
set Number of fields = 4
set Names of fields = upper, lower, mantle, seed
end
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = if(y>=80.e3, 1, 0); \
if(y<80.e3 && y>=70.e3, 1, 0); \
if(y<70.e3 && y>-100.e3,1, 0); \
if(y<68.e3 && y>60.e3 && x>=198.e3 && x<=202.e3 , 1, 0);
end
end
```

The initial thermal structure, radiogenic heating model and associated thermal properties are consistent with the prescribed thermal boundary conditions and produce a geotherm characteristic of the continental lithosphere. The equations defining the initial geotherm ([Chapman, 1986]) follow the form

where \(T\) is temperature, \(z\) is depth, \(T_T\) is the temperature at the layer surface (top), \(q_T\) is surface heat flux, \(k\) is thermal conductivity, and \(A\) is radiogenic heat production.

For a layer thickness \(\Delta z\), the basal temperature (\(T_B\)) and heat flux (\(q_B\))

In this example, specifying the top (273 \(K\)) temperature, surface heat flow (55 \(mW / m^2\)), thermal conductivity, and radiogenic heat production of each layer provides enough constraints to successively solve for the temperature and heat flux at the top of the lower crust and mantle.

As noted above, the initial zone of randomized plastic strain within the model center and strain softening of the friction and cohesion produces an initial pattern of distributed and unlocalized deformation across the zone of initial plastic strain (Fig. 102). After 5 million years of extension, distributed faulting is clearly evident in both the active and finite deformation fields and surface topography over an approximately 100 km wide region (Fig. 103). While deformation is distributed within this region, the fault system is clearly asymmetric as the rate of deformation and accumulated brittle strain varies between faults. Localization onto the two conjugate rift-bounding border faults is evident from the active deformation field. Notably, deformation of the free surface near the fixed left and right walls is evident at 5 Myr. Continued distortion of the surface mesh near the lateral boundaries may lead to solver issues, which can be overcome by either widening the model or allowing the mesh to deform along these boundaries.

With further extension for millions of years, significant crustal thinning and
surface topography development should occur in response to displacement along
the rift-bounding faults. However, given that the model only extends to 100 km
depth, the simulation will not produce a realistic representation of
continental breakup due to the lack of an upwelling asthenosphere layer.
Indeed, numerical studies that examine continental breakup, rather than just
the initial stages of continental extension, include an asthenospheric layer
or modified basal boundary conditions (e.g. Winkler boundary condition in [Brune *et al.*, 2014] for example) as temperature variations associated with
lithospheric thinning exert a first-order influence on the deformation
patterns. As noted below, numerous additional parameters may also affect the
temporal evolution of deformation patterns.