Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[WIP] Cookbook initial lithosphere composition, plastic strain, temperature and topography #5683

Open
wants to merge 34 commits into
base: main
Choose a base branch
from
Open

[WIP] Cookbook initial lithosphere composition, plastic strain, temperature and topography #5683

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
34 commits
Select commit Hold shift + click to select a range
35f9464
Initial plugins
anne-glerum Dec 12, 2023
46eda58
First stab prm file
anne-glerum Dec 12, 2023
3b477e5
Fix IC manager
anne-glerum Dec 12, 2023
26cf312
Set plastic strain grid to 10 km extra
anne-glerum Dec 12, 2023
1cfb54f
Move get_pointer to where it exists
anne-glerum Dec 13, 2023
aa833ba
Update prm
anne-glerum Dec 14, 2023
d45b26b
Store pointer and update read-in density initial temperature
anne-glerum Sep 30, 2022
a9ad749
Update read-in densities initial topography
anne-glerum Sep 30, 2022
a8cd5da
Remove unused parameter
anne-glerum Dec 17, 2023
6a95cae
Fix compiler warnings
anne-glerum Dec 17, 2023
088315c
Clean up
anne-glerum Jan 2, 2024
94f22a3
Clean up IT
anne-glerum Jan 2, 2024
0e769ec
Clean up IT
anne-glerum Jan 15, 2024
0cdf221
Fix type
anne-glerum Jan 15, 2024
cf8aaf8
Astyle
anne-glerum Jan 15, 2024
2ba7824
Indent
anne-glerum Jan 16, 2024
19ca0dd
Use function instead of cast
anne-glerum Jan 17, 2024
fabe1f0
Update comments topo
anne-glerum Jan 18, 2024
ae74ca5
Update topo files
anne-glerum Jan 19, 2024
2369d95
Fix rift-craton interaction
anne-glerum May 30, 2024
2fb1ceb
Use the IC function for local thicknesses
anne-glerum May 30, 2024
785dbdd
Fix overrides and parameter names
anne-glerum May 30, 2024
56c80ab
Update prm
anne-glerum May 30, 2024
eb89c81
Fix IC and IT
anne-glerum May 30, 2024
6847dd1
Rename prm
anne-glerum May 30, 2024
50b2a3e
Update text with figure
anne-glerum May 30, 2024
027213d
Add craton and AMR
anne-glerum May 30, 2024
5f26ef5
Update topo plugins
anne-glerum May 30, 2024
a051df2
Indent
anne-glerum May 30, 2024
dd54055
Add cookbook to geophysical setups
anne-glerum May 30, 2024
7b44330
Fix typos and md syntax
anne-glerum May 31, 2024
ae690a3
Add prm snippets and fix figure name
anne-glerum May 31, 2024
8e8239a
Expand description, add boundary conditions
anne-glerum Jun 17, 2024
bd49ff6
Add FS and traction BC
anne-glerum Jun 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
338 changes: 338 additions & 0 deletions cookbooks/initial_conditions_continental_rift/5p_fixed_CERI_craton_right.prm
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,338 @@
set Additional shared libraries = /Applications/ASPECT/VisualStudioCode/aspect/cookbooks/initial_conditions_continental_rift/libinitial_conditions_continental_rift.debug.so


set Dimension = 2
set End time = 3500
set Use years in output instead of seconds = true
set CFL number = 1
set Maximum time step = 3.5e3
set Maximum first time step = 2e2
set Maximum relative increase in time step = 100

set Adiabatic surface temperature = 1547
set Pressure normalization = no
set Output directory = /Applications/ASPECT/VisualStudioCode/aspect/cookbooks/initial_conditions_continental_rift/output

set Nonlinear solver scheme = single Advection, iterated Stokes
set Max nonlinear iterations = 2
set Nonlinear solver tolerance = 1e-3
set Max nonlinear iterations in pre-refinement = 0

# Solver parameters
subsection Solver parameters
subsection Stokes solver parameters
set Number of cheap Stokes solver steps = 1000
set Linear solver tolerance = 1e-3
set GMRES solver restart length = 50
set Use full A block as preconditioner = true
end
end

subsection Geometry model
set Model name = box

subsection Box
set X extent = 700e3
set Y extent = 300e3
set X repetitions = 70
set Y repetitions = 30
end

end

# Take the minimum of the continental geotherm
# and the adiabat with 'Adiabatic surface temperature'
subsection Initial temperature model
set List of model names = adiabatic, lithosphere with rift
set List of model operators = add, replace if valid
subsection Lithosphere with rift
set LAB isotherm temperature = 1623.
set Surface temperature = 293.
# Whether or not to prescribe the LAB depth up to a
# certain depth (in area below LAB)
# or use adiabat everywhere
set Use temperature compensation depth = true
set Temperature compensation depth = 200e3
end
subsection Adiabatic
# A reference profile of the compositional fields
# where x represents depth
subsection Function
# plastic_strain, viscous_strain, sediment_age, deposition_depth, sediment, upper, lower, mantle_L, ratio_marine_continental_sediment, silt_fraction
set Function expression = 0; \
0; \
0; \
0; \
0; \
if(x<23e3,1,0); \
if(x>=23e3&x<43e3,1,0); \
if(x>=43e3&x<200e3,1,0); \
0; \
0
end
set Age top boundary layer = 0
set Age bottom boundary layer = 0
end
end

subsection Compositional fields
set Number of fields = 10
set Names of fields = plastic_strain, viscous_strain, sediment_age, deposition_depth, sediment, upper, lower, mantle_L, ratio_marine_continental_sediment, silt_fraction
end

subsection Initial composition model
set List of model names = lithosphere with rift, rift box initial plastic strain, function
#set List of model names = lithosphere with rift, rift, function
subsection Lithosphere with rift
set Blend polygons and rifts = true
set Standard deviation of Gaussian rift geometry = 60e3
set Amplitude of Gaussian rift geometry = -0.25, 0, 0.17647
# LAB at 110 km
set Rift axis line segments = 350e3
set Layer thicknesses = 20e3, 15e3, 85e3

# Polygon LAB at 200 km
set Lithospheric polygon layer thicknesses = 23e3, 20e3 , 157.e3
set Lithospheric polygons = 450e3 > 3000e3
set Half width of polygon smoothing = 10e3
end

# Inherited noise
subsection Rift box initial plastic strain
set Standard deviation of Gaussian noise amplitude distribution = 200e3
set Maximum amplitude of Gaussian noise amplitude distribution = 0.250000
set Random number generator seed = 2
set Depth around which Gaussian noise is smoothed out = 50e3
set Halfwidth with which Gaussian noise is smoothed out in depth = 10e3
#set Grid intervals for noise X or radius = 2240
# Add 16 intervals for the additional 5 km that is added to the grid
# to account for initial topography
#set Grid intervals for noise Y or longitude = 976 #960
set Rift axis line segments = 350e3
end

subsection Function
set Variable names = x,y
set Function expression = 0; \
0; \
0; \
0; \
0; \
0; \
0; \
0; \
0; \
0
end
end

subsection Boundary composition model
set List of model names = function
set Fixed composition boundary indicators = bottom, top
subsection Function
set Variable names = x,y,t
set Function constants = h=300e3,My=1e-6, SL=0
# Names of fields = plastic_strain, viscous_strain, sediment_age, deposition_depth, sediment, upper, lower, mantle_L, ratio, silt_fraction
set Function expression = 0; \
0; \
if(y>h/2 && t>0,t*My,0); \
if(y>h/2 && t>0,h+SL-y,0); \
if(y>h/2 && t>0,1,0); \
if(y>h/2 && t==0,1,0); \
0; \
0; \
0; \
0
end
end

subsection Boundary temperature model
set List of model names = initial temperature, box
set List of model operators = add, minimum
set Fixed temperature boundary indicators = bottom, top
subsection Box
set Top temperature = 293
set Bottom temperature = 5000 # Unrealistically high, so that it is always taken from initial temperature plugin
end
end


subsection Mesh deformation
set Additional tangential mesh velocity boundary indicators = right
set Mesh deformation boundary indicators = top: free surface, top: ascii data

subsection Ascii data model
set Data directory = /Applications/ASPECT/VisualStudioCode/aspect/cookbooks/initial_conditions_continental_rift
set Data file name = 5p_fixed_craton_UC23km.txt
end
end

subsection Formulation
set Formulation = custom
# incompressible
set Mass conservation = ask material model
# use reference density in the temperature equation
set Temperature equation = reference density profile
end

subsection Boundary velocity model
set Prescribed velocity boundary indicators = left: function, right: function
set Tangential velocity boundary indicators =
subsection Function
set Coordinate system = cartesian
set Variable names = x,y,t
set Function constants = mm=0.001, yr=1, Z=300e3, X=700e3, outflow=10
set Function expression = if(x<X/2,-1,1)*outflow/2*mm/yr; \
0
end
end

subsection Boundary traction model
set Prescribed traction boundary indicators = bottom: initial lithostatic pressure
subsection Initial lithostatic pressure
set Representative point = 10000,0
end
end


subsection Gravity model
set Model name = vertical

subsection Vertical
set Magnitude = 9.81
end
end


subsection Material model
set Model name = visco plastic
set Material averaging = harmonic average only viscosity

subsection Visco Plastic

set Reference temperature = 293

set Minimum strain rate = 1.e-25
set Reference strain rate = 1.e-16

set Minimum viscosity = 5e18
set Maximum viscosity = 1e25

# bg, plastic_strain, viscous_strain, sediment_age, deposition_depth, sediment, upper, lower, mantle_L, ratio_marine_continental_sediment, silt_fraction
set Heat capacities = 1200
# Use densities for masked fields that are closest to the density of the area they are in to avoid artefacts from under/overshoot.
set Densities = 3182, 2613, 2758, 2410, 2410, 2410, 2613, 2758, 3163, 2410, 2410
set Thermal expansivities = 3.0e-5, 2.7e-5, 2.7e-5, 3.7e-5, 3.7e-5, 3.7e-5, 2.7e-5, 2.7e-5, 3.0e-5, 3.7e-5, 3.7e-5
set Define thermal conductivities = true
set Thermal conductivities = 3, 2.5, 2.5, 2.1, 2.1, 2.1, 2.5, 2.5, 3, 2.5, 2.5

set Viscosity averaging scheme = harmonic

set Viscous flow law = composite

# Dislocation creep parameters for
# 1. Background material/mantle (dry olivine)
# Hirth & Kohlstedt (2003), Geophys. Monogr. Am. Geophys. Soc., v.138, p.83-105.
# "Rheology of the upper mantle and the mantle wedge:a view from the experimentalists"
# 2. Upper crust (wet quartzite)
# Rutter & Brodie (2004), J. Struct. Geol., v.26, p.2011-2023.
# "Experimental grain size-sensitive flow of hot-pressed Brazilian quartz aggregates"
# 3. Lower crust and weak seed (wet anorthite)
# Rybacki et al. (2006), J. Geophys. Res., v.111(B3).
# "Influence of water fugacity and activation volume on the flow properties of fine-grained
# anorthite aggregates"
# Note that the viscous pre-factors below are scaled to plane strain from unixial strain experiments.
set Prefactors for dislocation creep = 2.12e-15, 1e-20, 1e-20 , 1e-20, 1e-20 , 8.57e-28, 8.57e-28, 7.13e-18, 6.52e-16, 1e-20, 1e-20
set Stress exponents for dislocation creep = 3.5, 1.0 , 1.0 , 1.0 , 1.0 , 4.0 , 4.0 , 3.0 , 3.5, 1, 1
set Activation energies for dislocation creep = 480.e3, 530.e3, 530.e3, 530.e3, 530.e3, 223.e3 , 223.e3 , 345.e3 , 530.e3, 530.e3, 530e3
set Activation volumes for dislocation creep = 11.e-6, 18.e-6, 18.e-6, 18.e-6, 18.e-6, 0. , 0. , 38.e-6 , 18.e-6, 18.e-6, 18e-6

set Prefactors for diffusion creep = 1.5e-9, 2.25e-9, 2.25e-9, 2.25e-9, 2.25e-9, 5.97e-19, 5.97e-19, 2.99e-25, 2.25e-9, 2.25e-9, 2.25e-9
set Activation energies for diffusion creep = 335.e3 , 375.e3 , 375.e3 , 375.e3 , 375.e3 , 223.e3 , 223.e3 , 159.e3 , 375.e3, 375e3, 375e3
set Activation volumes for diffusion creep = 4.e-6 , 6.e-6 , 6.e-6 , 6.e-6 , 6.e-6 , 0. , 0. , 38.e-6 , 6.e-6, 6.e-6, 6e-6
set Grain size = 1e-3
set Grain size exponents for diffusion creep = 0. , 0. , 0. , 0. , 0. , 2 , 2 , 3 , 0, 0, 0
# Plasticity parameters
set Angles of internal friction = 26.56
set Cohesions = 5.e6

# Strain weakening parameters
set Strain weakening mechanism = plastic weakening with plastic strain and viscous weakening with viscous strain
set Start plasticity strain weakening intervals = 0.0
set End plasticity strain weakening intervals = 2
set Cohesion strain weakening factors = 1.0
set Friction strain weakening factors = 0.25

set Start prefactor strain weakening intervals = 0.0
set End prefactor strain weakening intervals = 2
# TODO apply weakening to bg as well?
set Prefactor strain weakening factors = 1.0, 1.0, 1.0, 1.0, 1.0, 0.25, 0.25, 0.25, 0.25, 1.0, 1.0

end
end

subsection Mesh refinement
set Initial global refinement = 0
set Initial adaptive refinement = 0
set Time steps between mesh refinement = 0
set Strategy = minimum refinement function
set Skip solvers on initial refinement = false

subsection Minimum refinement function
set Coordinate system = cartesian
set Variable names = x,y
# Level 5 (highest) in upper 45 km (crust is at most 43 km thick).
# Level 3 up to 130 km depth (normal lithosphere is 120 km thick).
# Level 3 up to 210 km depth if craton present (craton is 200 km thick).
# Level 1 below.
set Function expression = if(y>=255e3,5,if(y>90e3, 3, 1))
end

end

subsection Heating model
set List of model names = compositional heating, adiabatic heating, shear heating
subsection Adiabatic heating
set Use simplified adiabatic heating = true
end
subsection Compositional heating
# order: background, plastic_strain, viscous_strain, sediment_age, deposition_depth, sediment, upper, lower, mantle_L, ratio, silt_fraction
set Use compositional field for heat production averaging = 0 , 0 , 0 , 0 , 0 , 1 , 1 , 1 , 1, 0, 0
set Compositional heating values = 0 , 0 , 0 , 0 , 0 , 1.2e-6 ,1.0e-6,0.1e-6, 0, 0, 0
end
end

subsection Adiabatic conditions model
set Model name = compute profile
subsection Compute profile
set Composition reference profile = function
# In terms of depth
# Moho depth craton 43 km, LAB depth 200 km (i.e. max LAB depth in model)
set Function expression = 0; \
0; \
0; \
0; \
0; \
if(x<23e3,1,0); \
if(x>=23e3&x<43e3,1,0); \
if(x>=43e3&x<200e3,1,0); \
0; \
0
end
end

subsection Postprocess
set List of postprocessors = velocity statistics, temperature statistics, visualization, composition statistics, topography
subsection Visualization
set Interpolate output = false
set Write higher order output = false
set List of output variables = adiabat, viscosity, density, named additional outputs, thermal conductivity, heating
set Time between graphical output = 500e3
set Point-wise stress and strain = true
end
subsection Topography
set Output to file = true
set Time between text output = 500e3
end

end
44 changes: 44 additions & 0 deletions cookbooks/initial_conditions_continental_rift/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,44 @@
# Copyright (C) 2015 - 2021 by the authors of the ASPECT code.
#
# This file is part of ASPECT.
#
# ASPECT is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2, or (at your option)
# any later version.
#
# ASPECT is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with ASPECT; see the file LICENSE. If not see
# <http://www.gnu.org/licenses/>.

CMAKE_MINIMUM_REQUIRED(VERSION 3.1.0)

FIND_PACKAGE(Aspect 2.4.0 QUIET HINTS ${Aspect_DIR} ../ ../../ $ENV{ASPECT_DIR})

IF (NOT Aspect_FOUND)
MESSAGE(FATAL_ERROR "\n"
"Could not find a valid ASPECT build/installation directory. "
"Please specify the directory where you are building ASPECT by passing\n"
" -D Aspect_DIR=<path to ASPECT>\n"
"to cmake or by setting the environment variable ASPECT_DIR in your shell "
"before calling cmake. See the section 'How to write a plugin' in the "
"manual for more information.")
ENDIF ()

DEAL_II_INITIALIZE_CACHED_VARIABLES()

SET(TARGET "initial_conditions_continental_rift")
PROJECT(${TARGET})

ADD_LIBRARY(${TARGET} SHARED
lithosphere_rift_IC.cc
lithosphere_rift_IT_2.cc
lithosphere_rift_topo.cc
rift_box_initial_plastic_strain.cc
)
ASPECT_SETUP_PLUGIN(${TARGET})
Loading
Loading