Modeling rifting#
For rifts to form, we need to localize deformation. In numerical models, this is often done using a seed that represents a weaker rheology compared to the surrounding material (e.g., a lower viscosity or a lower yield strength) or using a thermal anomaly.
We describe the process of rifting using two 2D models with viscoplastic rheology that are pulled apart along the side boundaries: one with a conductive initial temperature distribution and another with a thermal perturbation added to it (see figure below).
Initial strain rates in the model setups we use here to initiate rifting. Temperature profiles for the respective model is plotted on the right and the profile is marked by a white line on the left. Notice the increased strain rates in the model with an upwelling due to increased temperatures at those locations.
Geometry Model#
We use a rectangular box of 200x100 km
subsection Geometry model
set Model name = box
subsection Box
set X extent = 200000
set Y extent = 100000
end
end
Initial Compositional Field#
The compositions represent upper crust (20 km thick), lower crust (20 km thick) and mantle (60 km thick) are continuous horizontal layers of constant thickness. The non initial plastic strain is set to 0, while the initial plastic strain is randomized between 0.5 and 1.5. The randomized initial plastic strain physically represents a damage zone and is used to localize deformation in the model. Specifically, we reduce the cohesion and internal angle of frictional values in locations where the initial plastic strain value is between 0.5 and 1.5 by a factor of 0.25. The relevant weakening parameters are described in the Material model section.
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = 0; \
if(x>50.e3 && x<150.e3 && y>50.e3, 0.5 + rand_seed(1), 0); \
if(y>=80.e3, 1, 0); \
if(y<80.e3 && y>=60.e3, 1, 0); \
if(y<60.e3, 1, 0);
end
end
Location of the four compositional fields described in our model setup with background colors depicting the plastic_strain compositional value.
Initial Temperature Model#
The temperature field is represented by a linear increase with depth, with different thermal gradients representing different radiogenic heat production within each compositional layer. We have three layers: upper crust, lower crust and mantle, having different geothermal gradient within each layer (ith layer) determined by the radiogenic heat production (Ai), thermal conductivity (ki), and temperature (tsi) and heat flow (qsi) at top of that layer. The equation used to compute the temperature distribution is from eq. 4 from Chapman 1986.
Layer (i) |
|---|
upper crust (ts1, A1, k1, qs1) |
lower crust (ts2, A2, k2, qs2) |
mantle (ts3, A3, k3, qs3) |
subsection Initial temperature model
set List of model names = function
subsection Function
set Variable names = x,y
set Function constants = h=100e3, x0=200e3, ts1=273, ts2=633, ts3=893, \
A1=1.e-6, A2=0.25e-6, A3=0.0, r=20e3, \
k1=2.5, k2=2.5, k3=2.5, T_excess = 1800, \
qs1=0.055, qs2=0.035, qs3=0.030
set Function expression = if( (h-y)<=20.e3, \
ts1 + (qs1/k1)*(h-y) - (A1*(h-y)*(h-y))/(2.0*k1), \
if( (h-y)>20.e3 && (h-y)<=40.e3, \
ts2 + (qs2/k2)*(h-y-20.e3) - (A2*(h-y-20.e3)*(h-y-20.e3))/(2.0*k2), \
ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3) ) )
end
To the above temperatures, we add a gaussian thermal perturbation at the base in rifting_with_upwelling.prm setup to investigate the impact of upwelling on the rift formation as: amplitude*exp(-((x-x0/2)*(x-x0/2)+(y*y))/(2*d*d)), where amplitude is the excess temperature set to 500K, d is the width of the anomaly set to 20 km in our .prm file.
Material Model#
We use viscoplastic material model with different material properties for the upper crust, lower crust and mantle. Additionally, we weaken the cohesion and frictional coefficient values by a factor of 4 where the initial plastic strain is greater than 0.5 to promote localization and rifting in the model. The relevant parameter values are given below.
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
set Prefactors for dislocation creep = background: 7.37e-15, \
crust_upper: 1.37e-26, \
crust_lower: 5.71e-23, \
mantle_lithosphere: 7.37e-15
set Stress exponents for dislocation creep = background: 3.5, \
crust_upper: 4.0, \
crust_lower: 3.0, \
mantle_lithosphere: 3.5
set Activation energies for dislocation creep = background: 530.e3, \
crust_upper: 223.e3, \
crust_lower: 345.e3, \
mantle_lithosphere: 530.e3
set Activation volumes for dislocation creep = background: 18.e-6, \
crust_upper: 0, \
crust_lower: 0, \
mantle_lithosphere: 18.e-6
set Angles of internal friction = 30
set Cohesions = 20.e6
set Strain weakening mechanism = plastic weakening with plastic strain only
set Start plasticity strain weakening intervals = 0.5
set End plasticity strain weakening intervals = 1.5
set Cohesion strain weakening factors = 0.25
set Friction strain weakening factors = 0.25
end
end
We don’t reduce the cohesion and internal angle of friction in the model setup with a thermal upwelling, i.e., Cohesion strain weakening factors and Friction strain weakening factors are set to 1.
Boundary velocity model#
The side boundaries are pulled apart with a velocity of 2.5 mm/yr, and the equivalent inflow velocity is applied at the bottom boundary to maintain mass balance in the model. The top boundary is free surface.
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
