Hi,
Thanks a lot for the help John! Bit of a late response but I’ve spent a lot of time working on and off on these sorts of models. Here, I’ll show a simple 2D force driven continental extension model and then what I want to do with my project.
2D Force Driven Continental Extension Model
Full .prms and ascii data files - 2d_traction_1.prm (14.2 KB),
2d_traction_2.prm (14.3 KB),
continentalextension_2d_left.0.txt (20.5 KB),
continentalextension_2d_right.0.txt (20.9 KB)
As John suggested I used the initial lithostatic pressure condition with an additional traction function on the left and right boundaries to create a force that initiates rifting, I also applied a velocity y = 0 function so the model is only pulled in the x orientation, however, this isn’t necessary if you apply the traction force over a smaller area. The full boundary conditions I used were:
subsection Boundary velocity model
set Prescribed velocity boundary indicators = left y:function, right y:function
set Zero velocity boundary indicators = bottom
subsection Function
set Variable names = x,y
set Function expression = 0; 0
end
end
subsection Boundary traction model
set Prescribed traction boundary indicators = left:initial lithostatic pressure, right:initial lithostatic pressure, left:function, right:function
subsection Initial lithostatic pressure
set Representative point = 0e3, 0e3
end
subsection Function
set Variable names = x,y
set Function constants = stress= 2.7e8, point_1= 560e3, point_2= 520e3
set Function expression = if (x<100e3, if (y>340e3 && y<390e3, -stress, 0), if (y>340e3 && y<390e3, stress, 0)); 0
end
end
By applying a force like this the model produces high strain along the left and right boundaries (pictured below) where there are viscosity contrasts in the crust (i.e. across the upper-lower crust boundary). I’ve been avoiding this issue by only visualising the centre of the model, but if anyone has any advice on applying a force that won’t lead to these high strain effects let me know.
Furthermore, when the model starts necking, causing velocities to increase, I checkpoint the model and output the velocities along the left and right boundaries, then save them in an ascii data format. I then resume the model using this ascii data as set velocities (this is why there’s two .prm files).
I also spent a considerable amount of time working on getting the model to converge, the parameters I tweaked were:
set Nonlinear solver scheme = single Advection, iterated Stokes
set Nonlinear solver tolerance = 1e-4
set Max nonlinear iterations = 200
set CFL number = 0.3
set Maximum time step = 20e3
set Material averaging = harmonic average only viscosity
set Minimum viscosity = 1e18
set Maximum viscosity = 1e24
subsection Solver parameters
subsection Stokes solver parameters
set Stokes solver type = block AMG
set Number of cheap Stokes solver steps = 500
set Use full A block as preconditioner = true
set GMRES solver restart length = 100
end
subsection Newton solver parameters
set Maximum linear Stokes solver tolerance = 1e-2
set Use Eisenstat Walker method for Picard iterations = true
end
end
Early on the stokes solver converges in around 25 to 80 iterations, but after 100,000 years it only takes 4 to 10 iterations until necking where the number of iterations increase (full .log file:
log.txt). If anyone has any advice on how I could potentially improve the convergence of this model let me know.
My Project
For my project I want to use a 3D force driven model to determine how transform faults, velocities, and velocity vectors develop during changes in the traction orientation. To do this I start with a model that uses velocities on the left and right boundaries and extend the model at 1 cm/yr on each side until breakup (17 Myr). After breakup I then switch to using traction boundaries and change the traction vector to be oblique to the previous velocity angle. Additionally, the front and back boundaries are set to periodic.
Here, I tested applying tractions at various angles to see how the model’s velocities and velocity orientations would react after 17 Myr of extension using velocities:
A negative vector means pulling towards the top left, bottom right.
However, I’m a little confused by some of these results. Firstly, when I apply a traction vector at 0 degrees it produces a velocity vector of -32 degrees (which looks strange at first but is actually parallel to the transform) with a really low velocity of 0.1 cm/yr, but I’m a little unsure as to why the velocities are so low here and why they’re higher elsewhere, can anyone explain why the velocities are higher pulling obliquely at a negative angle? Additionally, there’s quite a large discrepancy between the traction and velocity vectors, with the velocity vectors being a lot larger, why is this? I presumed this is due to the free passage of material through the periodic boundaries and how the models already broken up, but the model’s quite sensitive to this effect and I wondered if anyone had any ideas to potentially reduce it? If you want any more information or have anything to add to my model setup let me know.
Thanks,
Luke
) - 
. I’m currently reading through some the new changes that have been made so I’m really excited to play around with it.
