diff --git a/cookbooks/initial_conditions_continental_rift/5p_fixed_CERI_craton_right.prm b/cookbooks/initial_conditions_continental_rift/5p_fixed_CERI_craton_right.prm new file mode 100644 index 00000000000..103eb58ea4c --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/5p_fixed_CERI_craton_right.prm @@ -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=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 diff --git a/cookbooks/initial_conditions_continental_rift/CMakeLists.txt b/cookbooks/initial_conditions_continental_rift/CMakeLists.txt new file mode 100644 index 00000000000..9d718c97571 --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/CMakeLists.txt @@ -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 +# . + +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=\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}) diff --git a/cookbooks/initial_conditions_continental_rift/continental_extension.prm b/cookbooks/initial_conditions_continental_rift/continental_extension.prm new file mode 100644 index 00000000000..fba92e9d2da --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/continental_extension.prm @@ -0,0 +1,310 @@ +#### Continental Extension Cookbook +# This cookbook is based off numerous published studies, five of which are listed below. +# For additional information, see these publications and references therein. +# 1. Brune, S., Heine, C., Perez-Gussinye, M., and Sobolev, S.V. (2014), Nat. Comm., v.5, n.4014, +# Rift migration explains continental margin asymmetry and crustal hyperextension +# 2. Huismans, R., and Beaumont, C. (2011), Nature, v.473, p.71-75. +# Depth-dependent extension, two-stage breakup and cratonic underplating at rifted margins +# 3. Naliboff, J., and Buiter, S.H. (2015), Earth Planet. Sci. Lett., v.421, p.58-67, +# "Rift Reactivation and migration during multiphase extension" +# 4. Naliboff, J., Glerum, A., Sascha, S., Peron-Pinvidic, G., and Wrona, T. (2020), Geophys. +# Res. Lett., 47, e2019GL086611, "Development of 3‐D rift heterogeneity through fault +# network evolution" +# 5. Sandiford, D., Brune, S., Glerum, A., Naliboff, J., and Whittaker, J.M. (2021), Geophys. +# Geochem. Geosys., 22, e2021GC009681, "Kinematics of footwall exhumation at oceanic +# detachment faults: Solid-block rotation and apparent unbending" + +#### Global parameters +set Additional shared libraries = /Applications/ASPECT/VisualStudioCode/aspect/cookbooks/initial_conditions_continental_rift/libinitial_conditions_continental_rift.debug.so +set Dimension = 2 +set Start time = 0 +set End time = 2e3 +set Use years in output instead of seconds = true +set Nonlinear solver scheme = single Advection, iterated Stokes +set Nonlinear solver tolerance = 1e-4 +set Max nonlinear iterations = 1 +set CFL number = 0.5 +set Maximum time step = 2e3 +set Output directory = output-continental_extension +set Pressure normalization = no + +#### Parameters describing the model + +# Governing equations +subsection Formulation + set Formulation = Boussinesq approximation +end + +# Model geometry (200x100 km, 20 km spacing) +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 70 + set Y repetitions = 30 + set X extent = 700e3 + set Y extent = 300e3 + end +end + +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Initial global refinement = 0 + 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 + +subsection Solver parameters + subsection Stokes solver parameters + set Number of cheap Stokes solver steps = 0 + set Linear solver tolerance = 1e-5 + end + +end + +# Advecting the free surface using a normal, rather than vertical, +# projection. To reduce mesh instabilities and associated solver +# issues when deformation becomes large, diffusion is applied to +# the free surface at each time step. +subsection Mesh deformation + set Mesh deformation boundary indicators = top: free surface + + subsection Free surface + set Surface velocity projection = normal + end + +end + +# Velocity on boundaries characterized by functions +# The outward velocity (x-direction) on the left and right walls is 0.25 cm/year +# The vertical velocity at the base is 0.25 cm/year (balances outflow on sides) +# Velocity components parallel to the base (x-velocity) and side walls (y-velocity) +# are unconstrained (i.e. 'free'). +subsection Boundary velocity model + set Prescribed velocity boundary indicators = left x: function, right x:function + + subsection Function + set Variable names = x,y + set Function constants = v=0.0025, w=700.e3, d=300.e3 + set Function expression = if (x < w/2 , -v, v) ; v*2*d/w + 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 +# Number and names of compositional fields +# The five compositional fields represent: +# 1. The plastic strain that accumualates over time, with the initial plastic strain removed +# 2. The plastic strain that accumulated over time, including the initial plastic strain values +# 3. The upper crust +# 4. The lower crust +# 5. The mantle lithosphere +subsection Compositional fields + set Number of fields = 5 + set Names of fields = noninitial_plastic_strain, plastic_strain, upper, lower, mantle_L +end + +# Initial values of different compositional fields +# The 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 and the initial plastic strain is randomized between +# 0.5 and 1.5. +subsection Initial composition model + set List of model names = rift box initial plastic strain, function + + 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.000000 + set Depth around which Gaussian noise is smoothed out = 50e3 + set Halfwidth with which Gaussian noise is smoothed out in depth = 10e3 + set Rift axis line segments = 350e3 + end + subsection Function + set Variable names = x,y + set Function expression = 0; \ + 0; \ + if(y>=280.e3, 1, 0); \ + if(y<280.e3 && y>=260.e3, 1, 0); \ + if(y<260.e3, 1, 0); + end +end + +# Composition: fixed on bottom (inflow boundary), free on sides and top +subsection Boundary composition model + set Fixed composition boundary indicators = bottom + set List of model names = initial composition +end + +# Temperature boundary conditions +# Top and bottom (fixed) temperatures are consistent with the initial temperature field +# Note that while temperatures are specified for the model sides, these values are +# not used as the sides are not specified "Fixed temperature boundaries". Rather, +# these boundaries are insulating (zero net heat flux). +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 + +# Initial temperature field +# Typical continental geotherm based on equations 4-6 from: +# D.S. Chapman (1986), "Thermal gradients in the continental crust", +# Geological Society of London Special Publications, v.24, p.63-70. +# The initial constraints are: +# Surface Temperature - upper crust (ts1) = 273 K +# Surface Heat Flow - upper crust (qs1) = 0.055 mW/m^2 +# Heat Production - upper crust (A1) = 1.00e-6 W/m^3; +# Heat Production - lower crust (A2) = 0.25e-6 W/m^3; +# Heat Production - mantle (A3) = 0.00e-6 W/m^3; +# Thermal Conductivity - all layers = 2.5 (W/(m K)); +# To satisfy these constraints, the following values are required: +# Surface Temperature - lower crust (ts2) = 633 K +# - mantle (ts3) = 893 K +# Surface Heat Flow - lower crust (qs2) = 0.035 W/m^2; +# - mantle (qs3) = 0.030 W/m^2; +# Note: The continental geotherm initial temperature model +# plugin can be used to compute an identical geotherm +# for the lithosphere. An example of how to use this +# plugin is illustrated in the test for this cookbook +# (tests/continental_extension.prm). +subsection Initial temperature model + set Model name = function + + subsection Function + set Variable names = x,y + set Function constants = h=100e3, ts1=273, ts2=633, ts3=893, \ + A1=1.e-6, A2=0.25e-6, A3=0.0, \ + k1=2.5, k2=2.5, k3=2.5, \ + 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 +end + +# Constant internal heat production values (W/m^3) for background material +# and compositional fields. +subsection Heating model + set List of model names = compositional heating + + subsection Compositional heating + set Use compositional field for heat production averaging = 1, 0, 0, 1, 1, 1 + set Compositional heating values = 0.0, 0.0, 0.0, 1.0e-6, 0.25e-6, 0. + end +end + +# Material model +# Rheology: Non-linear viscous flow and Drucker Prager Plasticity +# Values for most rheological parameters are specified for a background material and +# each compositional field. Values for viscous deformation are based on dislocation +# creep flow-laws, with distinct values for the upper crust (wet quartzite), lower +# crust (wet anorthite) and mantle (dry olivine). Table 1 of Naliboff and Buiter (2015), +# Earth Planet. Sci. Lett., v.421, p. 58-67 contains values for each of these flow laws. +subsection Material model + set Model name = visco plastic + + subsection Visco Plastic + # Reference temperature and viscosity + set Reference temperature = 273 + + # The minimum strain-rate helps limit large viscosities values that arise + # as the strain-rate approaches zero. + # The reference strain-rate is used on the first non-linear iteration + # of the first time step when the velocity has not been determined yet. + set Minimum strain rate = 1.e-20 + set Reference strain rate = 1.e-16 + + # Limit the viscosity with minimum and maximum values + set Minimum viscosity = 1e18 + set Maximum viscosity = 1e26 + + # Thermal diffusivity is adjusted to match thermal conductivities + # assumed in assigning the initial geotherm + set Define thermal conductivities = true + set Thermal conductivities = 2.5 + set Heat capacities = 750. + + # Density values of 1 are assigned to "strain" fields, which are not taken into + # account when computing material properties. + set Densities = 3300, 1.0, 1.0, 2700, 2900, 3300 + set Thermal expansivities = 2e-5 + + # Harmonic viscosity averaging + set Viscosity averaging scheme = harmonic + + # Choose to have the viscosity (pre-yield) follow a dislocation + # diffusion or composite flow law. Here, dislocation is selected + # so no need to specify diffusion creep parameters below, which are + # only used if "diffusion" or "composite" option is selected. + set Viscous flow law = dislocation + + # Dislocation creep parameters for + # 1. Background material/mantle (dry olivine) + # Hirth & Kohlstedt (2004), 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. + # For ease of identification, fields tracking strain are assigned prefactors of 1e-50 + set Prefactors for dislocation creep = 6.52e-16, 1.00e-50, 1.00e-50, 8.57e-28, 7.13e-18, 6.52e-16 + set Stress exponents for dislocation creep = 3.5, 1.0, 1.0, 4.0, 3.0, 3.5 + set Activation energies for dislocation creep = 530.e3, 0.0, 0.0, 223.e3, 345.e3, 530.e3 + set Activation volumes for dislocation creep = 18.e-6, 0.0, 0.0, 0.0, 0.0, 18.e-6 + + # Plasticity parameters + set Angles of internal friction = 30 + set Cohesions = 20.e6 + + # The parameters below weaken the friction and cohesion by a + # a factor of 4 between plastic strain values of 0.5 and 1.5. + 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 + +# Gravity model +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 9.81 + end +end + +# Post processing +subsection Postprocess + set List of postprocessors = composition statistics, heat flux densities, heat flux statistics, mass flux statistics, pressure statistics, temperature statistics, topography, velocity statistics, visualization + + subsection Visualization + set List of output variables = density, named additional outputs, strain rate, viscosity + set Output format = vtu + set Time between graphical output = 100.e3 + set Interpolate output = true + end +end diff --git a/cookbooks/initial_conditions_continental_rift/doc/Initial_density_with_craton.png b/cookbooks/initial_conditions_continental_rift/doc/Initial_density_with_craton.png new file mode 100644 index 00000000000..15b29aff3a6 Binary files /dev/null and b/cookbooks/initial_conditions_continental_rift/doc/Initial_density_with_craton.png differ diff --git a/cookbooks/initial_conditions_continental_rift/doc/Initial_plastic_strain_with_craton.png b/cookbooks/initial_conditions_continental_rift/doc/Initial_plastic_strain_with_craton.png new file mode 100644 index 00000000000..639b497c933 Binary files /dev/null and b/cookbooks/initial_conditions_continental_rift/doc/Initial_plastic_strain_with_craton.png differ diff --git a/cookbooks/initial_conditions_continental_rift/doc/Initial_temperature_with_craton.png b/cookbooks/initial_conditions_continental_rift/doc/Initial_temperature_with_craton.png new file mode 100644 index 00000000000..9e69ad3c51f Binary files /dev/null and b/cookbooks/initial_conditions_continental_rift/doc/Initial_temperature_with_craton.png differ diff --git a/cookbooks/initial_conditions_continental_rift/doc/Initial_topography_with_craton.png b/cookbooks/initial_conditions_continental_rift/doc/Initial_topography_with_craton.png new file mode 100644 index 00000000000..dd83c799b78 Binary files /dev/null and b/cookbooks/initial_conditions_continental_rift/doc/Initial_topography_with_craton.png differ diff --git a/cookbooks/initial_conditions_continental_rift/doc/initial_composition_lithosphere.prm b/cookbooks/initial_conditions_continental_rift/doc/initial_composition_lithosphere.prm new file mode 100644 index 00000000000..a34ad95ed99 --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/doc/initial_composition_lithosphere.prm @@ -0,0 +1,59 @@ +# Number and names of compositional fields +# The five compositional fields represent: +# 1. The plastic strain that accumualates over time, with the initial plastic strain removed +# 2. The plastic strain that accumulated over time, including the initial plastic strain values +# 3. The upper crust +# 4. The lower crust +# 5. The mantle lithosphere +subsection Compositional fields + set Number of fields = 5 + set Names of fields = noninitial_plastic_strain, plastic_strain, upper, lower, mantle_L +end + +# Initial values of different compositional fields +# The upper crust (20 km thick), lower crust (15 km thick) +# and mantle (85 km thick) are continuous horizontal layers +# of constant thickness, except in a center area. +# In the centre of the domain, the upper crust is thickened to 25 km, +# while the mantle lithosphere is thinned. The transition from this +# perturbed center to the reference thicknesses of the lithospheric +# layers is smoothened using a Gaussian. The combined effect +# is a hotter, weaker centre area that localizes deformation. +# An additional thicker/stronger area is added on the right of the domain +# to represent cratonic lithosphere. The transition in layer thicknesses +# uses a hyperbolic tangent. +# The non initial plastic strain is set +# to 0 and the initial plastic strain is randomized between +# 0.5 and 1.5 in an area around the rift axis. +subsection Initial composition model + set List of model names = lithosphere with rift, rift box initial plastic strain, function + + subsection Lithosphere with rift + # The reference layer thicknesses of the upper crust, lower crust and mantle lithosphere. + set Layer thicknesses = 20e3, 15e3, 85e3 + # How wide the weak zone will be + set Standard deviation of Gaussian rift geometry = 50e3 + # The maximum thinning/thickening of the lithospheric layers + set Amplitude of Gaussian rift geometry = -0.25, 0, 0.17647 + # Where the rift axis should be. In this case in the centre of the domain. + # In 3D simulations, line segments consisting of two coordinates + # should be provided. + set Rift axis line segments = 200e3 + + # The thicknesses of the lithospheric layer within a polygon + set Lithospheric polygon layer thicknesses = 20e3, 15e3 , 115.e3 + # The segments making up the polygon + set Lithospheric polygons = 550e3 > 3000e3 + # The half width of the transition from polygon to reference lithosphere + set Half width of polygon smoothing = 10e3 + end + + subsection Function + set Variable names = x,y + set Function expression = 0; \ + 0; \ + 0; \ + 0; \ + 0 + end +end diff --git a/cookbooks/initial_conditions_continental_rift/doc/initial_composition_strain.prm b/cookbooks/initial_conditions_continental_rift/doc/initial_composition_strain.prm new file mode 100644 index 00000000000..e56717798f2 --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/doc/initial_composition_strain.prm @@ -0,0 +1,22 @@ +subsection Initial composition model + set List of model names = lithosphere with rift, rift box initial plastic strain, function + + # Inherited plastic strain + subsection Rift box initial plastic strain + # We focus initial strain around the rift axis + set Rift axis line segments = 200e3 + # The strain values follow a Gaussian distribution around the rift axis + # The maximum magnitude of the strain + set Maximum amplitude of Gaussian noise amplitude distribution = 0.250000 + # The sigma width of the area of strain + set Standard deviation of Gaussian noise amplitude distribution = 50e3 + # The number with which to start of the random number generator + set Random number generator seed = 2 + # Using a hyperbolic tangent centred around this value, + # we smooth out the strain with increasing depth. + set Depth around which Gaussian noise is smoothed out = 50e3 + # The halfwidth of the hyperbolic tangent. + set Halfwidth with which Gaussian noise is smoothed out in depth = 10e3 + end + +end diff --git a/cookbooks/initial_conditions_continental_rift/doc/initial_conditions_continental_rift.md b/cookbooks/initial_conditions_continental_rift/doc/initial_conditions_continental_rift.md new file mode 100644 index 00000000000..70b45591b9e --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/doc/initial_conditions_continental_rift.md @@ -0,0 +1,158 @@ +# Initial and boundary conditions for continental rifting + +*This section was contributed by Anne Glerum.* + +Continental rifting is a process that thins the continental +lithosphere and can eventually lead to full break-up of that lithosphere +and seafloor spreading at a mid-ocean ridge {cite}[e.g., ][]`brune:2023`. +We can simulate continental rifting under far-field extensional +forces by prescribing certain velocities or stresses on the lateral +boundaries of our model domain (as in, for example, section +{ref}`sec:cookbooks:continental_extension`). +Here we choose to prescribe outward horizontal velocities on +the lateral boundaries. +In order to localize deformation under this prescribed +extension (and avoid localizing at the model boundaries), +we have to introduce initial weaknesses to the +model setup. The plugins prescribed in this cookbook set up +a lithospheric (upper crust, lower crust and mantle lithosphere) and asthenosphere +composition and temperature structure and then +add weaknesses by 1) smoothly perturbing the lithosphere +at user-specified locations and by 2) adding initial plastic +strain around these same locations. This strain represents +inherited weakness from for example prior deformation phases. +To represent other deviations of the reference lithosphere, +such as cratonic lithosphere, +additional regions of different lithospheric layer thicknesses +can be specified. +In case the perturbations of the lithosphere in combination with +a free surface or another type of mesh deformation induce large initial +velocities, an initial topography plugin approximates the topography resulting from +isostatic equilibrium with respect to the reference (unperturbed) lithosphere. + +Apart from prescribing tensional velocities on the lateral boundaries of +the box domain, we use the surface processes software FastScape to +modify the top boundary and a traction boundary condition for the bottom +boundary. This allows for the inclusion of sediment erosion, transport +and deposition, and the automatic balance of mass within the model domain. + +## Initial lithosphere composition and temperature + +```{figure-md} fig:setup + + +Lithospheric structure with perturbations that represent an +incipient rift and a cratonic region. +``` +To set up the reference lithospheric structure, we specify the Layer thicknesses +of the upper crust, lower crust and mantle lithosphere within the Lithosphere with rift +subsection (see input-file snippet below). The perturbations of these lithospheric layers - representing an +incipient rift and localizing deformation - follow Gaussian curves. +We can either thicken or thin the layers, depending on +the sign of each Gaussian amplitude specified: a negative amplitude thickens +the corresponding layer and vice versa. The width of the perturbation +is determined by the sigma of the Gaussians. The centre of the Gaussian +curve will be the rift centre and is set with the parameter Rift axis line segments: + + +```{literalinclude} initial_composition_lithosphere.prm +``` + +Note that you can specify multiple segments, both in 2D and in 3D. By changing +the thickness of the layers, we simultaneously change the temperature profile, +and thus the strength distribution, creating weaker lithosphere. + +In addition to an incipient rift, other deviations from the reference +lithospheric thicknesses can also be prescribed (e.g., representing thick, cold cratonic +lithosphere). In the above snippet, +these are specified as polygons (multiple polygons are allowed) with +their own layer thicknesses. The transition from reference lithosphere +to polygon lithosphere is implemented as a hyperbolic tangent with a +centre point defined by the coordinates of the polygon and a half-width (10 km in this case). + + +```{literalinclude} initial_temperature.prm +``` + + +For the initial temperature we combine two plugins: one that prescribes +a continental steady-state geotherm (lithosphere with rift) and one that computes an adiabatic +temperature profile (see snippet above). The geotherm implementation is based on Chapter 4.6 +of Turcotte and Schubert. For each x-coordinate, it solves a steady-state +heat conduction equation taking into account the respective densities, +heat productivities and thermal conductivities of the three lithospheric +layers and the user-set temperature at the top (Surface temperature) and +bottom of the lithosphere (LAB isotherm temperature). Below the lithosphere a NaN is returned, such that we can replace +this NaN by an adiabatic temperature profile. Care has to be taken to match +the user-specified LAB temperature and the adiabatic temperature at the same +depth by modifying the adiabatic surface temperature. One can choose to prescribe +the LAB temperature up to a certain compensation depth (e.g., the bottom of the +craton), such that the LAB and adiabatic temperatures match everywhere (Fig. {numref}`fig:initial_temperature`). + +```{figure-md} fig:initial_temperature + + + The initial temperature distribution. +``` +## Inherited heterogeneity + +```{literalinclude} initial_composition_strain.prm +``` + +To include inherited heterogeneity we also include initial plastic strain +around the rift centre (Fig. {numref}`fig:initial_strain`). This requires a +compositional field called plastic_strain and the plugin rift box initial plastic strain. +We generate random strain values between 0 and 1 that +are subsequently multiplied with a Gaussian distribution for which we specify +the maximum amplitude and sigma. To control the depth extent of the initial +strain, we specify a depth around which the amplitude is smoothed out to zero +with a hyperbolic tangent. The initial strain thus produced will be the same +for each run with the same specifications if the same number of MPI processes is used. +The plastic strain can be used to weaken the cohesion and/or internal angle of friction +of the plastic rheology, thus again providing weaknesses that localize deformation. + +```{figure-md} fig:initial_strain + + + The initial plastic strain distribution. +``` + +```{literalinclude} initial_topography.prm +``` + +## Initial topography + +Sometimes the combination of a free surface with lateral variations in the +lithospheric thicknesses leads to very high velocities in the first timesteps. +To avoid a drunken sailor effect and/or very small timesteps, we can then +prescribe initial topography that roughly fulfills isostatic equilibrium (Fig. {numref}`fig:initial_topography`). +The required topography is computed using the reference density for each material, +so it neglects temperature effects on the density. + + +```{figure-md} fig:initial_topography + + + The initial topography that roughly isostatically balances the rift perturbation and craton region + with respect to the reference lithosphere. +``` + +## Boundary conditions + +```{literalinclude} velocity_boundary_conditions.prm +``` +Extension is driven by outward horizontal velocities on the left and +right boundary of the model domain. To compensate for this outflow, +we open the bottom boundary to flow by prescribing a traction. The +prescribed traction is computed as the initial lithostatic pressure +at a user-specified point. We pick this point to lie in the reference +(unperturbed) lithosphere close to the left domain boundary. + +For the top boundary, we pick the mesh deformation plugin fastscape, +which provides an interface to the FastScape landscape evolution model. +This plugin computes changes to the mesh's surface topography based on +river incision, hill slope diffusion, sediment deposition and marine diffusion. +For all the different parameters relating to FastScape, have a look at +section {ref}`sec:cookbooks:fastscape_eroding_box`. We allow the top corner +of the right boundary to move with the surface, but not the left top corner point, +as this is our reference location. diff --git a/cookbooks/initial_conditions_continental_rift/doc/initial_temperature.prm b/cookbooks/initial_conditions_continental_rift/doc/initial_temperature.prm new file mode 100644 index 00000000000..26b74f5a00d --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/doc/initial_temperature.prm @@ -0,0 +1,50 @@ +# Temperature boundary conditions +# Top and bottom (fixed) temperatures are consistent with the initial temperature field. +# Because we use a free surface, the height/depth of the surface can change over time, +# meaning that if the initial temperature plugin is asked for the temperature +# at the surface, this temperature could also vary over time. Therefore +# we also use the box plugin and take the minimum of the initial temperature and +# the box plugins, thus maintaining a temperature of 293 K. +subsection Boundary temperature model + set Fixed temperature boundary indicators = bottom, top + set List of model names = initial temperature, box + set List of model operators = add, minimum + subsection Box + set Top temperature = 293 + # Unrealistically high, so that it is always taken from initial temperature plugin + set Bottom temperature = 5000 + end +end + +# Initial temperature field +# The initial temperature consists of a typical continental geotherm +# based on equations Turcotte and Schubert Ch. 4.6. +# Below the lithosphere-asthenosphere boundary (LAB, defined initially by the +# bottom of the mantle lithsphere and the 1623 K isotherm, +# we prescribe an adiabatic temperature increase. +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 = 1573. + set Surface temperature = 293. + # Whether or not to prescribe the LAB temperature up to a + # certain depth (in the area below LAB) + # or to use adiabat everywhere below the LAB. + set Use temperature compensation depth = true + set Temperature compensation depth = 150e3 + end + subsection Adiabatic + # A reference profile of the compositional fields + # where x represents depth + subsection Function + # noninitial_plastic_strain, plastic_strain, upper, lower, mantle_L + set Function expression = 0; \ + 0; \ + if(x<20e3,1,0); \ + if(x>=20e3&x<35e3,1,0); \ + if(x>=35e3&x<120e3,1,0) + end + end + +end diff --git a/cookbooks/initial_conditions_continental_rift/doc/initial_topography.prm b/cookbooks/initial_conditions_continental_rift/doc/initial_topography.prm new file mode 100644 index 00000000000..34ef795c2e9 --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/doc/initial_topography.prm @@ -0,0 +1,19 @@ +# Model geometry (400x200 km, 40 km spacing) +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 10 + set Y repetitions = 5 + set X extent = 400e3 + set Y extent = 200e3 + end + + # With respect to the reference lithosphere, + # move the surface of perturbed lithosphere + # to obtain isostasy. + subsection Initial topography model + set Model name = lithosphere with rift + end + +end diff --git a/cookbooks/initial_conditions_continental_rift/doc/velocity_boundary_conditions.prm b/cookbooks/initial_conditions_continental_rift/doc/velocity_boundary_conditions.prm new file mode 100644 index 00000000000..826efdc814f --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/doc/velocity_boundary_conditions.prm @@ -0,0 +1,27 @@ +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= 155e3, 4, 3) + end +end + + +subsection Solver parameters + subsection Stokes solver parameters + set Number of cheap Stokes solver steps = 0 + end +end + +# Use FastScape to deform the mesh. +subsection Mesh deformation + set Additional tangential mesh velocity boundary indicators = right + set Mesh deformation boundary indicators = top: fastscape + + subsection Free surface + set Free surface stabilization theta = 1 + end + + subsection Fastscape + set Vertical exaggeration = -1 + set Maximum timestep length = 500 + set Number of fastscape timesteps per aspect timestep = 10 + set Maximum surface refinement level = 4 + set Surface refinement difference = 0 + set Additional fastscape refinement = 0 + set Use ghost nodes = true + set Y extent in 2d = 25e3 + set Uplift and advect with fastscape = true + + subsection Boundary conditions + set Front = 1 + set Right = 1 + set Back = 1 + set Left = 1 + end + + subsection Erosional parameters + set Drainage area exponent = 0.4 #m + set Slope exponent = 1 #n + set Multi-direction slope exponent = -1 #p + + set Bedrock diffusivity = 5e-3 #kd + set Bedrock river incision rate = 1e-5 #kf + set Bedrock deposition coefficient = 1 #G + + set Sediment diffusivity = -1 + set Sediment river incision rate = -1 + set Sediment deposition coefficient = 1 #G + end + + set Use marine component = true + set Sediment rain rates = 1e-4 # m/yr + set Sediment rain time intervals = + + # No difference between sand and silt + subsection Marine parameters + set Sea level = -200 + set Sand porosity = 0. + set Silt porosity = 0. + set Sand e-folding depth = 3000 # m + set Silt e-folding depth = 3000 # m + + set Sand-silt ratio = 0.5 + set Depth averaging thickness = 1e3 # m + set Sand transport coefficient = 150 # m2/yr + set Silt transport coefficient = 150 # m2/yr + end + end +end + +# Velocity on boundaries characterized by functions +# The outward velocity (x-direction) on the left and right walls is 0.25 cm/year +# The vertical velocity at the base is 0.25 cm/year (balances outflow on sides) +# Velocity components parallel to the base (x-velocity) and side walls (y-velocity) +# are unconstrained (i.e. 'free'). +subsection Boundary velocity model + set Prescribed velocity boundary indicators = left: function, right: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) ; 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 + +# Number and names of compositional fields +# The five compositional fields represent: +# 1. The plastic strain that accumualates over time, with the initial plastic strain removed +# 2. The plastic strain that accumulated over time, including the initial plastic strain values +# 3. The upper crust +# 4. The lower crust +# 5. The mantle lithosphere +subsection Compositional fields + set Number of fields = 5 + set Names of fields = noninitial_plastic_strain, plastic_strain, upper, lower, mantle_L +end + +# Initial values of different compositional fields +# The upper crust (20 km thick), lower crust (15 km thick) +# and mantle (85 km thick) are continuous horizontal layers +# of constant thickness, except in a center area. +# In the centre of the domain, the upper crust is thickened to 25 km, +# while the mantle lithosphere is thinned. The transition from this +# perturbed center to the reference thicknesses of the lithospheric +# layers is smoothened using a Gaussian. The combined effect +# is a hotter, weaker centre area that localizes deformation. +# An additional thicker/stronger area is added on the right of the domain +# to represent cratonic lithosphere. The transition in layer thicknesses +# uses a hyperbolic tangent. +# The non initial plastic strain is set +# to 0 and the initial plastic strain is randomized between +# 0.5 and 1.5 in an area around the rift axis. +subsection Initial composition model + set List of model names = lithosphere with rift, rift box initial plastic strain, function + + subsection Lithosphere with rift + # The reference layer thicknesses of the upper crust, lower crust and mantle lithosphere. + set Layer thicknesses = 20e3, 15e3, 85e3 + # How wide the weak zone will be + set Standard deviation of Gaussian rift geometry = 50e3 + # The maximum thinning/thickening of the lithospheric layers + set Amplitude of Gaussian rift geometry = -0.25, 0, 0.17647 + # Where the rift axis should be. In this case in the centre of the domain. + # In 3D simulations, line segments consisting of two coordinates + # should be provided. + set Rift axis line segments = 200e3 + + # The thicknesses of the lithospheric layer within a polygon + set Lithospheric polygon layer thicknesses = 20e3, 15e3 , 115.e3 + # The segments making up the polygon + set Lithospheric polygons = 350e3 > 3000e3 + # The half width of the transition from polygon to reference lithosphere + set Half width of polygon smoothing = 10e3 + end + + # Inherited plastic strain + subsection Rift box initial plastic strain + # We focus initial strain around the rift axis + set Rift axis line segments = 200e3 + # The strain values follow a Gaussian distribution around the rift axis + # The maximum magnitude of the strain + set Maximum amplitude of Gaussian noise amplitude distribution = 0.250000 + # The sigma width of the area of strain + set Standard deviation of Gaussian noise amplitude distribution = 50e3 + # The number with which to start of the random number generator + set Random number generator seed = 2 + # Using a hyperbolic tangent centred around this value, + # we smooth out the strain with increasing depth. + set Depth around which Gaussian noise is smoothed out = 50e3 + # The halfwidth of the hyperbolic tangent. + set Halfwidth with which Gaussian noise is smoothed out in depth = 10e3 + end + + subsection Function + set Variable names = x,y + set Function expression = 0; \ + 0; \ + 0; \ + 0; \ + 0 + end +end + +# Composition: fixed on bottom (inflow boundary), free on sides and top +subsection Boundary composition model + set Fixed composition boundary indicators = bottom + set List of model names = function + subsection Function + set Function expression = 0; \ + 0; \ + 0; \ + 0; \ + 0 + end +end + +# Temperature boundary conditions +# Top and bottom (fixed) temperatures are consistent with the initial temperature field. +# Because we use a free surface, the height/depth of the surface can change over time, +# meaning that if the initial temperature plugin is asked for the temperature +# at the surface, this temperature could also vary over time. Therefore +# we also use the box plugin and take the minimum of the initial temperature and +# the box plugins, thus maintaining a temperature of 293 K. +subsection Boundary temperature model + set Fixed temperature boundary indicators = bottom, top + set List of model names = initial temperature, box + set List of model operators = add, minimum + subsection Box + set Top temperature = 293 + # Unrealistically high, so that it is always taken from initial temperature plugin + set Bottom temperature = 5000 + end +end + +# Initial temperature field +# The initial temperature consists of a typical continental geotherm +# based on equations Turcotte and Schubert Ch. 4.6. +# Below the lithosphere-asthenosphere boundary (LAB, defined initially by the +# bottom of the mantle lithsphere and the 1623 K isotherm, +# we prescribe an adiabatic temperature increase. +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 = 1573. + set Surface temperature = 293. + end + subsection Adiabatic + # A reference profile of the compositional fields + # where x represents depth + subsection Function + # noninitial_plastic_strain, plastic_strain, upper, lower, mantle_L + set Function expression = 0; \ + 0; \ + if(x<20e3,1,0); \ + if(x>=20e3&x<35e3,1,0); \ + if(x>=35e3&x<120e3,1,0) + end + end + +end + +# Constant internal heat production values (W/m^3) for background material +# and compositional fields. +subsection Heating model + set List of model names = compositional heating + + subsection Compositional heating + set Use compositional field for heat production averaging = 0, 0, 0, 1, 1, 1 + set Compositional heating values = 0.0, 0.0, 0.0, 1.0e-6, 0.25e-6, 0. + end +end + +# Material model +# Rheology: Non-linear viscous flow and Drucker Prager Plasticity +# Values for most rheological parameters are specified for a background material and +# each compositional field. Values for viscous deformation are based on dislocation +# creep flow-laws, with distinct values for the upper crust (wet quartzite), lower +# crust (wet anorthite) and mantle (dry olivine). Table 1 of Naliboff and Buiter (2015), +# Earth Planet. Sci. Lett., v.421, p. 58-67 contains values for each of these flow laws. +subsection Material model + set Model name = visco plastic + + subsection Visco Plastic + # Reference temperature and viscosity + set Reference temperature = 273 + + # The minimum strain-rate helps limit large viscosities values that arise + # as the strain-rate approaches zero. + # The reference strain-rate is used on the first non-linear iteration + # of the first time step when the velocity has not been determined yet. + set Minimum strain rate = 1.e-20 + set Reference strain rate = 1.e-16 + + # Limit the viscosity with minimum and maximum values + set Minimum viscosity = 1e18 + set Maximum viscosity = 1e26 + + # Thermal diffusivity is adjusted to match thermal conductivities + # assumed in assigning the initial geotherm + set Define thermal conductivities = true + set Thermal conductivities = 2.5 + set Heat capacities = 750. + + # Density values of 1 are assigned to "strain" fields, which are not taken into + # account when computing material properties. + set Densities = 3300, 1.0, 1.0, 2700, 2900, 3200 + set Thermal expansivities = 2e-5 + + # Harmonic viscosity averaging + set Viscosity averaging scheme = harmonic + + # Choose to have the viscosity (pre-yield) follow a dislocation + # diffusion or composite flow law. Here, dislocation is selected + # so no need to specify diffusion creep parameters below, which are + # only used if "diffusion" or "composite" option is selected. + set Viscous flow law = dislocation + + # Dislocation creep parameters for + # 1. Background material/mantle (dry olivine) + # Hirth & Kohlstedt (2004), 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. + # For ease of identification, fields tracking strain are assigned prefactors of 1e-50 + set Prefactors for dislocation creep = 6.52e-16, 1.00e-50, 1.00e-50, 8.57e-28, 7.13e-18, 6.52e-16 + set Stress exponents for dislocation creep = 3.5, 1.0, 1.0, 4.0, 3.0, 3.5 + set Activation energies for dislocation creep = 530.e3, 0.0, 0.0, 223.e3, 345.e3, 530.e3 + set Activation volumes for dislocation creep = 18.e-6, 0.0, 0.0, 0.0, 0.0, 18.e-6 + + # Plasticity parameters + set Angles of internal friction = 30 + set Cohesions = 20.e6 + + # The parameters below weaken the friction and cohesion by a + # a factor of 4 between plastic strain values of 0.5 and 1.5. + 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 + +# Gravity model +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 9.81 + end +end + +# Post processing +subsection Postprocess + set List of postprocessors = composition statistics, temperature statistics, topography, velocity statistics, visualization + + subsection Visualization + set List of output variables = density, named additional outputs, strain rate, viscosity + set Output format = vtu + set Time between graphical output = 100.e3 + set Interpolate output = true + end +end diff --git a/cookbooks/initial_conditions_continental_rift/initial_conditions_continental_rift_no_initial_topography.prm b/cookbooks/initial_conditions_continental_rift/initial_conditions_continental_rift_no_initial_topography.prm new file mode 100644 index 00000000000..33007a01a4f --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/initial_conditions_continental_rift_no_initial_topography.prm @@ -0,0 +1,385 @@ +#### This is an adaptation of the Continental Extension Cookbook +#### using new plugins for the initial composition, temperature +#### topography and strain of an incipient continental rift system. + +#### Global parameters +set Additional shared libraries = $ASPECT_SOURCE_DIR/cookbooks/initial_conditions_continental_rift/libinitial_conditions_continental_rift.debug.so +set Dimension = 2 +set Start time = 0 +set End time = 0 +set Use years in output instead of seconds = true +set Nonlinear solver scheme = single Advection, iterated Stokes +set Nonlinear solver tolerance = 1e-5 +set Max nonlinear iterations = 50 +set CFL number = 0.5 +set Maximum time step = 20e3 +set Output directory = output-initial_conditions_continental_rift_no_initial_topography_with_craton +set Pressure normalization = no +set Adiabatic surface temperature = 1573 + +#### Parameters describing the model + +# Governing equations +subsection Formulation + set Formulation = Boussinesq approximation +end + +# Model geometry (400x200 km, 40 km spacing) +subsection Geometry model + set Model name = box + + subsection Box + set X repetitions = 10 + set Y repetitions = 5 + set X extent = 400e3 + set Y extent = 200e3 + end +end + +# Globally refine the mesh to 5 km spacing, and then +# adaptively refine the mesh to 2.5 km spacing for the crustal layers. +# These values ensure areas +# undergoing brittle deformation are in the high-resolution +# region. +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>= 155e3, 4, 3) + end +end + + +subsection Solver parameters + subsection Stokes solver parameters + set Number of cheap Stokes solver steps = 0 + end +end + + +# Use FastScape to deform the mesh. +subsection Mesh deformation + set Additional tangential mesh velocity boundary indicators = right + set Mesh deformation boundary indicators = top: fastscape + + subsection Free surface + set Free surface stabilization theta = 1 + end + + subsection Fastscape + set Vertical exaggeration = -1 + set Maximum timestep length = 500 + set Number of fastscape timesteps per aspect timestep = 10 + set Maximum surface refinement level = 4 + set Surface refinement difference = 0 + set Additional fastscape refinement = 0 + set Use ghost nodes = true + set Y extent in 2d = 25e3 + set Uplift and advect with fastscape = true + + subsection Boundary conditions + set Front = 1 + set Right = 1 + set Back = 1 + set Left = 1 + end + + subsection Erosional parameters + set Drainage area exponent = 0.4 #m + set Slope exponent = 1 #n + set Multi-direction slope exponent = -1 #p + + set Bedrock diffusivity = 5e-3 #kd + set Bedrock river incision rate = 1e-5 #kf + set Bedrock deposition coefficient = 1 #G + + set Sediment diffusivity = -1 + set Sediment river incision rate = -1 + set Sediment deposition coefficient = 1 #G + end + + set Use marine component = true + set Sediment rain rates = 1e-4 # m/yr + set Sediment rain time intervals = + + # No difference between sand and silt + subsection Marine parameters + set Sea level = -200 + set Sand porosity = 0. + set Silt porosity = 0. + set Sand e-folding depth = 3000 # m + set Silt e-folding depth = 3000 # m + + set Sand-silt ratio = 0.5 + set Depth averaging thickness = 1e3 # m + set Sand transport coefficient = 150 # m2/yr + set Silt transport coefficient = 150 # m2/yr + end + end +end + +# Velocity on boundaries characterized by functions +# The outward velocity (x-direction) on the left and right walls is 0.25 cm/year +# The vertical velocity at the base is 0.25 cm/year (balances outflow on sides) +# Velocity components parallel to the base (x-velocity) and side walls (y-velocity) +# are unconstrained (i.e. 'free'). +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 + +# Number and names of compositional fields +# The five compositional fields represent: +# 1. The plastic strain that accumualates over time, with the initial plastic strain removed +# 2. The plastic strain that accumulated over time, including the initial plastic strain values +# 3. The upper crust +# 4. The lower crust +# 5. The mantle lithosphere +subsection Compositional fields + set Number of fields = 5 + set Names of fields = noninitial_plastic_strain, plastic_strain, upper, lower, mantle_L +end + +# Initial values of different compositional fields +# The upper crust (20 km thick), lower crust (15 km thick) +# and mantle (85 km thick) are continuous horizontal layers +# of constant thickness, except in a center area. +# In the centre of the domain, the upper crust is thickened to 25 km, +# while the mantle lithosphere is thinned. The transition from this +# perturbed center to the reference thicknesses of the lithospheric +# layers is smoothened using a Gaussian. The combined effect +# is a hotter, weaker centre area that localizes deformation. +# An additional thicker/stronger area is added on the right of the domain +# to represent cratonic lithosphere. The transition in layer thicknesses +# uses a hyperbolic tangent. +# The non initial plastic strain is set +# to 0 and the initial plastic strain is randomized between +# 0.5 and 1.5 in an area around the rift axis. +subsection Initial composition model + set List of model names = lithosphere with rift, rift box initial plastic strain, function + + subsection Lithosphere with rift + # The reference layer thicknesses of the upper crust, lower crust and mantle lithosphere. + set Layer thicknesses = 20e3, 15e3, 85e3 + # How wide the weak zone will be + set Standard deviation of Gaussian rift geometry = 50e3 + # The maximum thinning/thickening of the lithospheric layers + set Amplitude of Gaussian rift geometry = -0.25, 0, 0.17647 + # Where the rift axis should be. In this case in the centre of the domain. + # In 3D simulations, line segments consisting of two coordinates + # should be provided. + set Rift axis line segments = 200e3 + + # The thicknesses of the lithospheric layer within a polygon + set Lithospheric polygon layer thicknesses = 20e3, 15e3 , 115.e3 + # The segments making up the polygon + set Lithospheric polygons = 350e3 > 3000e3 + # The half width of the transition from polygon to reference lithosphere + set Half width of polygon smoothing = 10e3 + end + + # Inherited plastic strain + subsection Rift box initial plastic strain + # We focus initial strain around the rift axis + set Rift axis line segments = 200e3 + # The strain values follow a Gaussian distribution around the rift axis + # The maximum magnitude of the strain + set Maximum amplitude of Gaussian noise amplitude distribution = 0.250000 + # The sigma width of the area of strain + set Standard deviation of Gaussian noise amplitude distribution = 50e3 + # The number with which to start of the random number generator + set Random number generator seed = 2 + # Using a hyperbolic tangent centred around this value, + # we smooth out the strain with increasing depth. + set Depth around which Gaussian noise is smoothed out = 50e3 + # The halfwidth of the hyperbolic tangent. + set Halfwidth with which Gaussian noise is smoothed out in depth = 10e3 + end + + subsection Function + set Variable names = x,y + set Function expression = 0; \ + 0; \ + 0; \ + 0; \ + 0 + end +end + +# Composition: fixed on bottom (inflow boundary), free on sides and top +subsection Boundary composition model + set Fixed composition boundary indicators = bottom + set List of model names = function + subsection Function + set Function expression = 0; \ + 0; \ + 0; \ + 0; \ + 0 + end +end + +# Temperature boundary conditions +# Top and bottom (fixed) temperatures are consistent with the initial temperature field. +# Because we use a free surface, the height/depth of the surface can change over time, +# meaning that if the initial temperature plugin is asked for the temperature +# at the surface, this temperature could also vary over time. Therefore +# we also use the box plugin and take the minimum of the initial temperature and +# the box plugins, thus maintaining a temperature of 293 K. +subsection Boundary temperature model + set Fixed temperature boundary indicators = bottom, top + set List of model names = initial temperature, box + set List of model operators = add, minimum + subsection Box + set Top temperature = 293 + # Unrealistically high, so that it is always taken from initial temperature plugin + set Bottom temperature = 5000 + end +end + +# Initial temperature field +# The initial temperature consists of a typical continental geotherm +# based on equations Turcotte and Schubert Ch. 4.6. +# Below the lithosphere-asthenosphere boundary (LAB, defined initially by the +# bottom of the mantle lithsphere and the 1623 K isotherm, +# we prescribe an adiabatic temperature increase. +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 = 1573. + set Surface temperature = 293. + end + subsection Adiabatic + # A reference profile of the compositional fields + # where x represents depth + subsection Function + # noninitial_plastic_strain, plastic_strain, upper, lower, mantle_L + set Function expression = 0; \ + 0; \ + if(x<20e3,1,0); \ + if(x>=20e3&x<35e3,1,0); \ + if(x>=35e3&x<120e3,1,0) + end + end + +end + +# Constant internal heat production values (W/m^3) for background material +# and compositional fields. +subsection Heating model + set List of model names = compositional heating + + subsection Compositional heating + set Use compositional field for heat production averaging = 0, 0, 0, 1, 1, 1 + set Compositional heating values = 0.0, 0.0, 0.0, 1.0e-6, 0.25e-6, 0. + end +end + +# Material model +# Rheology: Non-linear viscous flow and Drucker Prager Plasticity +# Values for most rheological parameters are specified for a background material and +# each compositional field. Values for viscous deformation are based on dislocation +# creep flow-laws, with distinct values for the upper crust (wet quartzite), lower +# crust (wet anorthite) and mantle (dry olivine). Table 1 of Naliboff and Buiter (2015), +# Earth Planet. Sci. Lett., v.421, p. 58-67 contains values for each of these flow laws. +subsection Material model + set Model name = visco plastic + + subsection Visco Plastic + # Reference temperature and viscosity + set Reference temperature = 273 + + # The minimum strain-rate helps limit large viscosities values that arise + # as the strain-rate approaches zero. + # The reference strain-rate is used on the first non-linear iteration + # of the first time step when the velocity has not been determined yet. + set Minimum strain rate = 1.e-20 + set Reference strain rate = 1.e-16 + + # Limit the viscosity with minimum and maximum values + set Minimum viscosity = 1e18 + set Maximum viscosity = 1e26 + + # Thermal diffusivity is adjusted to match thermal conductivities + # assumed in assigning the initial geotherm + set Define thermal conductivities = true + set Thermal conductivities = 2.5 + set Heat capacities = 750. + + # Density values of 1 are assigned to "strain" fields, which are not taken into + # account when computing material properties. + set Densities = 3300, 1.0, 1.0, 2700, 2900, 3300 + set Thermal expansivities = 2e-5 + + # Harmonic viscosity averaging + set Viscosity averaging scheme = harmonic + + # Choose to have the viscosity (pre-yield) follow a dislocation + # diffusion or composite flow law. Here, dislocation is selected + # so no need to specify diffusion creep parameters below, which are + # only used if "diffusion" or "composite" option is selected. + set Viscous flow law = dislocation + + # Dislocation creep parameters for + # 1. Background material/mantle (dry olivine) + # Hirth & Kohlstedt (2004), 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. + # For ease of identification, fields tracking strain are assigned prefactors of 1e-50 + set Prefactors for dislocation creep = 6.52e-16, 1.00e-50, 1.00e-50, 8.57e-28, 7.13e-18, 6.52e-16 + set Stress exponents for dislocation creep = 3.5, 1.0, 1.0, 4.0, 3.0, 3.5 + set Activation energies for dislocation creep = 530.e3, 0.0, 0.0, 223.e3, 345.e3, 530.e3 + set Activation volumes for dislocation creep = 18.e-6, 0.0, 0.0, 0.0, 0.0, 18.e-6 + + # Plasticity parameters + set Angles of internal friction = 30 + set Cohesions = 20.e6 + + # The parameters below weaken the friction and cohesion by a + # a factor of 4 between plastic strain values of 0.5 and 1.5. + 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 + +# Gravity model +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 9.81 + end +end + +# Post processing +subsection Postprocess + set List of postprocessors = composition statistics, temperature statistics, topography, velocity statistics, visualization + + subsection Visualization + set List of output variables = density, named additional outputs, strain rate, viscosity + set Output format = vtu + set Time between graphical output = 100.e3 + set Interpolate output = true + end +end diff --git a/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IC.cc b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IC.cc new file mode 100644 index 00000000000..a929d933c4b --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IC.cc @@ -0,0 +1,404 @@ +/* + Copyright (C) 2017 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 + . +*/ + + +#include "lithosphere_rift_IC.h" +#include "lithosphere_rift_IT.h" +#include +#include + + +namespace aspect +{ + namespace InitialComposition + { + template + void + LithosphereRift:: + initialize () + { + // Check that the corresponding initial temperature model is used. + const std::vector active_initial_temperature_models = this->get_initial_temperature_manager().get_active_initial_temperature_names(); + AssertThrow(find(active_initial_temperature_models.begin(),active_initial_temperature_models.end(), "lithosphere with rift") != active_initial_temperature_models.end(), + ExcMessage("The lithosphere with rift initial composition plugin requires the lithosphere with rift initial temperature plugin.")); + } + + template + double + LithosphereRift:: + initial_composition (const Point &position, + const unsigned int compositional_index) const + { + // Retrieve the indices of the fields that represent the lithospheric layers. + // We assume a 3-layer system with an upper crust, lower crust and lithospheric mantle. + const unsigned int id_upper = this->introspection().compositional_index_for_name("upper"); + const unsigned int id_lower = this->introspection().compositional_index_for_name("lower"); + const unsigned int id_mantle_L = this->introspection().compositional_index_for_name("mantle_L"); + + // Determine coordinate system + const bool cartesian_geometry = Plugins::plugin_type_matches> (this->get_geometry_model()); + + // Get the surface coordinates of the point under consideration + const Point surface_point = surface_position(position, cartesian_geometry); + + // Compute the local thickness of the upper crust, lower crust and mantle part of the lithosphere + // (in this exact order) based on the distance from the rift axis and the polygons. + const std::vector local_thicknesses = compute_local_thicknesses(surface_point); + + // Get depth with respect to the surface. + const double depth = this->get_geometry_model().depth(position); + + // Check which layer the current point lies in and return a value of 1 if the field corresponds to the layer. + if (depth <= local_thicknesses[0] && compositional_index == id_upper) + return 1.; + else if (depth > local_thicknesses[0] && depth <= local_thicknesses[0] + local_thicknesses[1] + && compositional_index == id_lower) + return 1.; + else if (depth > local_thicknesses[0] + local_thicknesses[1] && depth <= local_thicknesses[0] + local_thicknesses[1] + local_thicknesses[2] + && compositional_index == id_mantle_L) + return 1.; + else + return 0.; + } + + template + double + LithosphereRift:: + distance_to_rift (const Point &surface_position) const + { + // Initiate distance with large value + double distance_to_rift_axis = 1e23; + double temp_distance = 0; + + // Loop over all line segments + for (unsigned int i_segments = 0; i_segments < rift_point_list.size(); ++i_segments) + { + // The Utilities function only works in 3d, so compute distance here in 2d. + temp_distance = (dim == 2) ? std::abs(surface_position[0]-rift_point_list[i_segments][0][0]) + : std::abs(Utilities::distance_to_line(rift_point_list[i_segments], Point<2>(surface_position[0],surface_position[dim-2]))); + + // Get the minimum distance + distance_to_rift_axis = std::min(distance_to_rift_axis, temp_distance); + } + + return distance_to_rift_axis; + } + + template + Point + LithosphereRift:: + surface_position (const Point &position, + const bool cartesian_geometry) const + { + Point surface_point; + if (cartesian_geometry) + { + for (unsigned int d=0; d spherical_point = Utilities::Coordinates::cartesian_to_spherical_coordinates(position); + // return lon [degrees], lat [degrees] + for (unsigned int d=0; d + std::pair + LithosphereRift:: + distance_to_polygon (const Point &surface_position) const + { + // Inside the polygon is positive, outside negative. We assume the different polygons do not overlap. + // The Utilities function only works in 3d, so compute distance here in 2d. + double max_distance = -1e24; + unsigned int max_distance_polygon = 0; + for (unsigned int n = 0; npolygon_point_list[n][0][0] && surface_position[0](polygon_point_list[n], Point<2>(surface_position[0],surface_position[dim-2])); + } + + if (temp_distance > max_distance) + { + max_distance = temp_distance; + max_distance_polygon = n; + } + } + return std::pair (max_distance, max_distance_polygon); + } + + template + std::vector + LithosphereRift:: + compute_local_thicknesses(const Point &surface_point) const + { + // Get the distance to the rift segments along a path parallel to the surface + const double distance_to_rift_axis = distance_to_rift(surface_point); + + // Get the signed distance to potential polygons of different lithospheric thicknesses + const std::pair distance_to_L_polygon = distance_to_polygon(surface_point); + + // Compute the local thickness of the upper crust, lower crust and mantle part of the lithosphere + // (in this exact order) based on the distance from the rift axis and the polygons. + // The transition from reference and/or rift perturbation lithosphere to polygon lithosphere + // is smoothed by a hyperbolic tangent. + const double polygon_contribution = (0.5 + 0.5 * std::tanh(distance_to_L_polygon.first / sigma_polygon)); + const double rift_contribution = (0.5 - 0.5 * std::tanh(distance_to_L_polygon.first / sigma_polygon)); + + // The rift perturbation follows a Gaussian contribution. + std::vector local_thicknesses(3); + for (unsigned int i = 0; i < 3; ++i) + local_thicknesses[i] = (1.0 - A_rift[i] * std::exp((-std::pow(distance_to_rift_axis, 2) / (2.0 * std::pow(sigma_rift, 2))))) + * reference_thicknesses[i] * rift_contribution + polygon_contribution * polygon_thicknesses[distance_to_L_polygon.second][i]; + + return local_thicknesses; + } + + template + void + LithosphereRift::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection ("Initial composition model"); + { + prm.enter_subsection("Lithosphere with rift"); + { + prm.declare_entry("Standard deviation of Gaussian rift geometry", "20000", + Patterns::Double(0), + "The standard deviation of the Gaussian distribution of the thinning/thickening " + "of the lithosphere thicknesses. This parameter is taken to be the same " + "for all rift segments. Units: \\si{\\meter} or degrees."); + prm.declare_entry("Amplitude of Gaussian rift geometry", "0.2", + Patterns::List(Patterns::Double(-1, 1)), + "The amplitude of the Gaussian distribution of the thinning/thickening of the. " + "lithosphere thicknesses. This parameter is taken to be the same for all rift segments, " + "but can be set to vary per lithosphere layer. " + "Units: none."); + prm.declare_entry("Half width of polygon smoothing", "20000", + Patterns::Double(0), + "The half width of the hyperbolic tangent smoothing used to transition to the " + "lithospheric thicknesses of the polygon. This parameter is taken to be the same for all polygons. " + "Units: \\si{\\meter} or degrees."); + prm.declare_entry("Layer thicknesses", "30000.", + Patterns::List(Patterns::Double(0)), + "List of reference lithospheric layer thicknesses, i.e., the thicknesses of " + "the upper crust, lower crust and lithospheric mantle layers. " + "Units: \\si{\\meter}"); + prm.declare_entry("Rift axis line segments", + "", + Patterns::Anything(), + "The line segments that represent the rift axis separated by a semi-colon. " + "Each segment is made up of two points that represent horizontal coordinates (x,y) or (lon,lat). " + "The exact format for the point list describing the segments is " + "\"x1,y1>x2,y2;x2,y2>x3,y3;x4,y4>x5,y5\". Note that the segments can be connected " + "or isolated. The units of the coordinates are " + "dependent on the geometry model. In the box model they are in meters, in the " + "chunks they are in degrees. Units: \\si{\\meter} or degrees."); + prm.declare_entry("Lithospheric polygons", + "", + Patterns::List(Patterns::Anything()), + "The points making up polygons that represent an area of different lithospheric thickness. " + "The polygons are separated by semicolons. Each polygon is a list of " + "points that represent horizontal coordinates (x,y) or (lon,lat). " + "The exact format for the point list describing a polygon is " + "\"x1,y1>x2,y2>x3,y3>x4,y4>x5,y5\". Note that the polygon is assumed to be closed. " + "The units of the coordinates are dependent on the geometry model. " + "In the box model they are in meters, in the chunks they are in degrees. Units: \\si{\\meter} or degrees."); + prm.declare_entry("Lithospheric polygon layer thicknesses", "30000.", + Patterns::List(Patterns::List(Patterns::Double(0), 0, 3, ","), 0, 10, ";"), + "List of thicknesses of the lithospheric layers for each polygon." + "For each polygon, a total of 3 thicknesses should be given (upper crust, lower crust, mantle lithosphere)." + "If only one value is given, then all layers are assigned the same value. Units: \\si{\\meter}"); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + template + void + LithosphereRift::parse_parameters (ParameterHandler &prm) + { + // Check that the required compositional fields exist. + AssertThrow(this->introspection().compositional_name_exists("upper"), + ExcMessage("We need a compositional field called 'upper' representing the upper crust.")); + AssertThrow(this->introspection().compositional_name_exists("lower"), + ExcMessage("We need a compositional field called 'lower' representing the lower crust.")); + AssertThrow(this->introspection().compositional_name_exists("mantle_L"), + ExcMessage("We need a compositional field called 'mantle_L' representing the lithospheric part of the mantle.")); + + prm.enter_subsection ("Initial composition model"); + { + prm.enter_subsection("Lithosphere with rift"); + { + sigma_rift = prm.get_double ("Standard deviation of Gaussian rift geometry"); + A_rift = Utilities::possibly_extend_from_1_to_N(Utilities::string_to_double(Utilities::split_string_list(prm.get("Amplitude of Gaussian rift geometry"))), + 3, + "Amplitude of Gaussian rift geometry"); + sigma_polygon = prm.get_double ("Half width of polygon smoothing"); + reference_thicknesses = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Layer thicknesses"))), + 3, + "Layer thicknesses"); + // Read in the string of segments + const std::string temp_all_segments = prm.get("Rift axis line segments"); + // Split the string into segment strings + const std::vector temp_segments = Utilities::split_string_list(temp_all_segments,';'); + const unsigned int n_temp_segments = temp_segments.size(); + rift_point_list.resize(n_temp_segments); + // Loop over the segments to extract the points + for (unsigned int i_segment = 0; i_segment < n_temp_segments; ++i_segment) + { + // In 3d a line segment consists of 2 points, + // in 2d of only 1 (the horizontal ridge axis runs orthogonal to x and y into the screen) + const std::vector temp_segment = Utilities::split_string_list(temp_segments[i_segment],'>'); + + if (dim == 3) + { + AssertThrow(temp_segment.size() == 2,ExcMessage ("The given segment '" + temp_segment[i_segment] + "' is not correct. " + "It should only contain 2 parts: " + "the two points of the segment, separated by a '>'.")); + } + else + { + AssertThrow(temp_segment.size() == 1,ExcMessage ("The given segment '" + temp_segment[i_segment] + "' is not correct. " + "It should only contain 1 part: " + "the point representing the rift axis.")); + } + + // Loop over the dim-1 points of each segment + for (unsigned int i_points = 0; i_points < dim-1; ++i_points) + { + const std::vector temp_point = Utilities::string_to_double(Utilities::split_string_list(temp_segment[i_points],',')); + + if (dim == 3) + { + AssertThrow(temp_point.size() == 2,ExcMessage ("The given coordinates of segment '" + temp_segment[i_points] + "' are not correct. " + "It should only contain 2 parts: " + "the two coordinates of the segment end point, separated by a ','.")); + } + else + { + AssertThrow(temp_point.size() == 1,ExcMessage ("The given coordinates of segment '" + temp_segment[i_points] + "' are not correct. " + "It should only contain 1 part: " + "the one coordinate of the segment end point.")); + } + + // Add the point to the list of points for this segment. Two points are set + // even in 2d, although in 2d they are the same. + rift_point_list[i_segment][i_points][0] = temp_point[0]; + rift_point_list[i_segment][i_points][1] = temp_point[dim-2]; + } + } + + // Split the polygon string into the separate polygons + const std::vector temp_polygons = Utilities::split_string_list(prm.get("Lithospheric polygons"),';'); + const std::vector temp_thicknesses = Utilities::split_string_list(prm.get("Lithospheric polygon layer thicknesses"),';'); + const unsigned int n_polygons = temp_polygons.size(); + AssertThrow(temp_thicknesses.size() == n_polygons, + ExcMessage("The number of polygons (" + Utilities::int_to_string(n_polygons) + + ") does not correspond to the number of polygons for which a thickness is prescribed (" + + Utilities::int_to_string(temp_thicknesses.size()) + ").")); + + polygon_point_list.resize(n_polygons); + polygon_thicknesses.resize(n_polygons); + + for (unsigned int i_polygons = 0; i_polygons < n_polygons; ++i_polygons) + { + polygon_thicknesses[i_polygons] = Utilities::string_to_double(Utilities::split_string_list(temp_thicknesses[i_polygons],',')); + AssertThrow(polygon_thicknesses[i_polygons].size() == 3, + ExcMessage ("The number of layer thicknesses should be equal to 3 for polygon: " + Utilities::int_to_string(i_polygons) + + " but it is " + Utilities::int_to_string(polygon_thicknesses[i_polygons].size()))); + + // Split the polygon string into point strings + const std::vector temp_points = Utilities::split_string_list(temp_polygons[i_polygons],'>'); + const unsigned int n_temp_points = temp_points.size(); + if (dim == 3) + { + AssertThrow(n_temp_points>=3, ExcMessage ("The number of polygon points should be equal to or larger than 3 in 3d.")); + } + else + { + AssertThrow(n_temp_points==2, ExcMessage ("The number of polygon points should be equal to 2 in 2d.")); + } + + polygon_point_list[i_polygons].resize(n_temp_points); + + // Loop over the points of the polygon. Each point coordinate should consist of 2 (in 3d) or 1 (in 2d) values. + for (unsigned int i_points = 0; i_points < n_temp_points; ++i_points) + { + const std::vector temp_point = Utilities::string_to_double(Utilities::split_string_list(temp_points[i_points],',')); + AssertThrow(temp_point.size() == dim-1,ExcMessage ("The given coordinates of point '" + temp_points[i_points] + "' are not correct. " + "It should only contain 1 (2d) or 2 (in 3d) parts: " + "the longitude/x (and latitude/y in 3d) coordinate (separated by a ',').")); + + // Add the point to the list of points for this segment + polygon_point_list[i_polygons][i_points][0] = temp_point[0]; + polygon_point_list[i_polygons][i_points][1] = temp_point[dim-2]; + } + + // For simplicity later on, we want the polygon coordinates of each point in 2d to be increasing. + if (dim == 2) + AssertThrow(polygon_point_list[i_polygons][0][0] < polygon_point_list[i_polygons][1][0], ExcMessage("The order of the x coordinates of the 2 points " + "of each 2d polygon should be ascending. ")); + + } + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + } +} + +// explicit instantiations +namespace aspect +{ + namespace InitialComposition + { + ASPECT_REGISTER_INITIAL_COMPOSITION_MODEL(LithosphereRift, + "lithosphere with rift", + "A class that implements initial conditions for a continental rift " + "by computing the variable thickness of the upper crust, lower crust and mantle " + "lithosphere layers. Thicknesses can change towards the axis of the " + "rift according to a Gaussian distribution whose amplitude and standard " + "deviation can be specified from the input file. " + "Additional variations in layer thicknesses can be introduced through the " + "addition of polygons. Note that the different polygons " + "are assumed not to overlap. The user can decide how overlapping rifts " + "and polygons are treated. ") + } +} diff --git a/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IC.h b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IC.h new file mode 100644 index 00000000000..80cfb48be7a --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IC.h @@ -0,0 +1,153 @@ +/* + Copyright (C) 2023 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 + . +*/ + + +#ifndef _aspect_initial_composition_lithosphere_rift_h +#define _aspect_initial_composition_lithosphere_rift_h + +#include +#include + +namespace aspect +{ + namespace InitialComposition + { + using namespace dealii; + + /** + * A class that implements initial conditions for the composition + * considering a lithosphere consisting of an upper crust, + * lower crust and mantle part. + * + * @ingroup InitialCompositionModels + */ + template + class LithosphereRift : public Interface, + public SimulatorAccess + { + public: + + /** + * Initialization function. This function is called once at the + * beginning of the program. Checks preconditions. + */ + void + initialize () override; + + /** + * Return the initial composition as a function of position and number + * of compositional field. + */ + virtual + double initial_composition (const Point &position, + const unsigned int compositional_index) const override; + + /** + * Return the overall shortest distance to the rift center segments. + */ + double distance_to_rift (const Point &position) const; + + /* + * Return the overall shortest distance to the polygon segments. Polygons + * can be used to specify lithosphere of different thicknesses. + */ + std::pair distance_to_polygon (const Point &position) const; + + /* + * Return the position of the point in surface coordinates, + * i.e. x(,y) in meters or lon(,lat) in degrees. + */ + Point surface_position (const Point &position, + const bool cartesian_geometry) const; + + /* + * Compute the local thicknesses of the upper and lower crust + * and lithospheric mantle based on the distance to the rift center + * and the edge of polygons. + */ + std::vector compute_local_thicknesses(const Point &position) const; + + /** + * Declare the parameters this class takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + virtual + void + parse_parameters (ParameterHandler &prm) override; + + private: + + /** + * The standard deviation of the rift-perpendicular Gaussian distribution + * of the thinning/thickening of the lithospheric thicknesses with its + * maximum along the rift axis. + */ + double sigma_rift; + + /** + * The maximum amplitude of the Gaussian distribution of the thinning/thickening + * of the lithospheric thicknesses with distance from the rift axis. + * The amplitude should have values between -1 and 1, where positive + * numbers represent a reduction in thickness and negative numbers an increase. + * For example, values of 0.25, 0, 0 reduce the reference thickness of the + * upper crust by 25%, while the lower crust and mantle lithosphere are + * untouched. + */ + std::vector A_rift; + + /** + * The half width of the hyperbolic tangent used to smooth the transitions + * between reference and polygon lithospheric thicknesses. + */ + double sigma_polygon; + + /** + * The list of line segments consisting of two 2d coordinates per segment. + * The segments represent the rift axis. + */ + std::vector,2 >> rift_point_list; + + /** + * The list of lists of polygon points. + * The polygons can represent areas of different lithospheric thicknesses. + */ + std::vector>> polygon_point_list; + + /** + * Vector for the reference field thicknesses away from rift and polygons. + */ + std::vector reference_thicknesses; + + /** + * Vector for the field thicknesses for each polygon. + */ + std::vector> polygon_thicknesses; + }; + } +} + + +#endif diff --git a/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IT.h b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IT.h new file mode 100644 index 00000000000..17d26e8af0b --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IT.h @@ -0,0 +1,188 @@ +/* + Copyright (C) 2014 - 2016 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 doc/COPYING. If not see + . +*/ + + + +#ifndef _aspect_initial_temperature_lithosphere_rift_h +#define _aspect_initial_temperature_lithosphere_rift_h + +#include +#include +#include + +namespace aspect +{ + + + namespace InitialTemperature + { + + /** + * + * @ingroup InitialTemperatures + */ + template + class LithosphereRift : public Interface, public ::aspect::SimulatorAccess + { + public: + /** + * Constructor. + */ + LithosphereRift (); + + /** + * Initialization function. This function is called once at the + * beginning of the program. Checks preconditions. + */ + void + initialize () override; + + /** + * Return the initial temperature as a function of position. + */ + virtual + double initial_temperature (const Point &position) const override; + + /** + * Return the initial temperature as a function of depth and + * the local layer thicknesses. + */ + virtual + double temperature (const double depth, + const std::vector thicknesses) const; + + /** + * Declare the parameters this class takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + virtual + void + parse_parameters (ParameterHandler &prm) override; + + private: + + /** + * Surface temperature + */ + double T0; + + /** + * LAB isotherm temperature + */ + double LAB_isotherm; + + /** + * Vector for lithospheric layer heat production rates. + */ + std::vector heat_productivities; + + /** + * Vector for lithospheric layer thermal conductivities. + */ + std::vector conductivities; + + /** + * Vector for lithospheric layer densities. + */ + std::vector densities; + + /** + * Vector for the reference field thicknesses. + */ + std::vector reference_thicknesses; + + /** + * The standard deviation of the Gaussian distribution of the lithospheric thicknesses + * around the rift segments. + */ + double sigma_rift; + + /** + * The maximum amplitude of the Gaussian distribution of the thinning/thickening + * of the lithospheric thicknesses with distance from the rift axis. + * The amplitude should have values between -1 and 1, where positive + * numbers represent a reduction in thickness and negative numbers an increase. + * For example, values of 0.25, 0, 0 reduce the reference thickness of the + * upper crust by 25%, while the lower crust and mantle lithosphere are + * untouched. + */ + std::vector A_rift; + + /** + * Vector for the polygon thicknesses. + */ + std::vector> polygon_thicknesses; + + /** + * The half width of the hyperbolic tangent used to smooth the transitions + * between reference and polygon lithospheric thicknesses. + */ + double sigma_polygon; + + /** + * Whether or not to take the polygon thicknesses as dominant, or to smooth them + * gradually into rift areas. + */ + bool blend_rift_and_polygon; + + /** + * Whether or not to use a compensation depth for the temperature + * up to which the LAB isotherm is prescribed. If yes, temperatures + * in the sub-lithospheric mantle are set to the LAB isotherm up + * to the depth held by the parameter compensation_depth. + */ + bool use_compensation_depth; + double compensation_depth; + + /** + * A shared pointer to the initial composition object + * that ensures that the current object can continue + * to access the initial composition object beyond the + * first time step. + */ + std::shared_ptr> initial_composition_manager; + + /** + * Vector containing the number of phases for each composition. + */ + std::unique_ptr> n_phases_for_each_composition; + + /** + * Vector containing the names of the compositional fields. + */ + std::vector list_of_composition_names; + + /** + * Whether a cartesian box geometry is used, or a geometry + * of spherical type (e.g., shell or chunk). + */ + bool cartesian_geometry; + }; + } +} + + +#endif diff --git a/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IT_2.cc b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IT_2.cc new file mode 100644 index 00000000000..38cf202c95d --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_IT_2.cc @@ -0,0 +1,316 @@ +/* + Copyright (C) 2011 - 2016 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 doc/COPYING. If not see + . +*/ + + +#include "lithosphere_rift_IT.h" +#include +#include "lithosphere_rift_IC.h" +#include +#include +#include +#include + +#include + +namespace aspect +{ + namespace InitialTemperature + { + template + LithosphereRift::LithosphereRift () + {} + + + template + void + LithosphereRift:: + initialize () + { + + // Check that the required radioactive heating model ("compositional heating") is used. + const std::vector &heating_models = this->get_heating_model_manager().get_active_heating_model_names(); + AssertThrow(std::find(heating_models.begin(), heating_models.end(), "compositional heating") != heating_models.end(), + ExcMessage("The lithosphere with rift initial temperature plugin requires the compositional heating plugin.")); + + // Check that the required material model ("visco plastic") is used. + AssertThrow(Plugins::plugin_type_matches> (this->get_material_model()), + ExcMessage("The lithosphere with rift initial temperature plugin requires the viscoplastic material model plugin.")); + + // The simulator only keeps the initial conditions around for + // the first time step. As a consequence, we have to save a + // shared pointer to that object ourselves the first time we get here. + if (initial_composition_manager == nullptr) + initial_composition_manager = this->get_initial_composition_manager_pointer(); + + // Check that the required initial composition model is used. + // This doesn't work during initialization. + //AssertThrow(initial_composition_manager->template has_matching_initial_composition_model>(), + // ExcMessage("The initial temperature plugin lithosphere with rift requires the corresponding initial composition plugin.")); + + // Determine whether a cartesian or a spherical geometry is used. + cartesian_geometry = Plugins::plugin_type_matches> (this->get_geometry_model()); + } + + + template + double + LithosphereRift:: + initial_temperature (const Point &position) const + { + Assert(initial_composition_manager->template has_matching_initial_composition_model>(), + ExcMessage("The initial temperature plugin lithosphere with rift requires the corresponding initial composition plugin.")); + + // Get the initial composition plugin + const InitialComposition::LithosphereRift &ic = initial_composition_manager->template get_matching_initial_composition_model>(); + + // Convert to surface position + const Point surface_position = ic.surface_position(position, cartesian_geometry); + + // Compute the local thickness of the upper crust, lower crust and mantle part of the lithosphere + // based on the distance from the rift axis and polygon edges. + std::vector local_thicknesses = ic.compute_local_thicknesses(surface_position); + + // Get depth with respect to the surface. + const double depth = this->get_geometry_model().depth(position); + + return temperature(depth, local_thicknesses); + } + + + template + double + LithosphereRift:: + temperature (const double depth, + const std::vector layer_thicknesses) const + { + // Compute some constants + const double a = 0.5*densities[0]*heat_productivities[0]*layer_thicknesses[0] + 0.5*densities[1]*heat_productivities[1]*layer_thicknesses[1] + conductivities[0]/layer_thicknesses[0]*T0; + const double b = 1./(conductivities[0]/layer_thicknesses[0]+conductivities[1]/layer_thicknesses[1]); + const double c = 0.5*densities[1]*heat_productivities[1]*layer_thicknesses[1] + conductivities[2]/layer_thicknesses[2]*LAB_isotherm; + const double d = 1./(conductivities[1]/layer_thicknesses[1]+conductivities[2]/layer_thicknesses[2]); + + //Temperature at boundary between layer 1 and 2 + const double T1 = (a*b + conductivities[1]/layer_thicknesses[1]*c*d*b) / (1.-(conductivities[1]*conductivities[1])/(layer_thicknesses[1]*layer_thicknesses[1])*d*b); + //Temperature at boundary between layer 2 and 3 + const double T2 = (c + conductivities[1]/layer_thicknesses[1]*T1) * d; + + // Temperature in layer 1 + if (depth < layer_thicknesses[0]) + return -0.5*densities[0]*heat_productivities[0]/conductivities[0]*std::pow(depth,2) + (0.5*densities[0]*heat_productivities[0]*layer_thicknesses[0]/conductivities[0] + (T1-T0)/layer_thicknesses[0])*depth + T0; + // Temperature in layer 2 + else if (depth < layer_thicknesses[0]+layer_thicknesses[1]) + return -0.5*densities[1]*heat_productivities[1]/conductivities[1]*std::pow(depth-layer_thicknesses[0],2.) + (0.5*densities[1]*heat_productivities[1]*layer_thicknesses[1]/conductivities[1] + (T2-T1)/layer_thicknesses[1])*(depth-layer_thicknesses[0]) + T1; + // Temperature in layer 3 + else if (depth < layer_thicknesses[0]+layer_thicknesses[1]+layer_thicknesses[2]) + return (LAB_isotherm-T2)/layer_thicknesses[2] *(depth-layer_thicknesses[0]-layer_thicknesses[1]) + T2; + // Temperature in the sublithospheric mantle up to the user-set compensation depth + // equals the temperature at the lithosphere-asthenosphere boundary. + else if (use_compensation_depth && depth < compensation_depth) + return LAB_isotherm; + // Return a NaN sublithospheric temperature. + // This way we can combine the continental geotherm with an adiabatic profile from the input file + // using the "if valid" operator on the "List of initial temperature models". + else + return std::numeric_limits::quiet_NaN(); + } + + + template + void + LithosphereRift::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection ("Initial temperature model"); + { + prm.enter_subsection("Lithosphere with rift"); + { + prm.declare_entry("Surface temperature", "273.15", + Patterns::Double(0), + "The value of the surface temperature. Units: \\si{\\kelvin}."); + prm.declare_entry("LAB isotherm temperature", "1673.15", + Patterns::Double(0), + "The value of the isothermal boundary temperature assumed at the LAB " + "and up to the compensation depth when 'Use temperature compensation depth' " + "is set to true. Units: \\si{\\kelvin}."); + prm.declare_entry ("Use temperature compensation depth", "false", + Patterns::Bool (), + "Whether or not to use a compensation depth up to which the LAB isotherm temperature " + "is prescribed below the lithosphere. If true, this plugin can be combined with " + "the adabiatic temperature plugin. The adiabatic surface temperature can be adapted to " + "ensure the LAB temperature and the adiabatic temperature match at the compensation depth. " + "If false, combining this plugin with the adiabatic temperature plugin will lead to " + "jumps in the temperature profile at the LAB when the lithospheric thickness varies " + "laterally. Units: -."); + prm.declare_entry("Temperature compensation depth", "120e3", + Patterns::Double(0), + "The depth of the temperature compensation, i.e. the depth up to which the LAB isotherm " + "is prescribed below the lithosphere. Units: \\si{\\meter}."); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + + + template + void + LithosphereRift::parse_parameters (ParameterHandler &prm) + { + unsigned int n_fields = 0; + prm.enter_subsection ("Compositional fields"); + { + n_fields = prm.get_integer ("Number of fields"); + } + prm.leave_subsection(); + + prm.enter_subsection ("Initial temperature model"); + { + prm.enter_subsection("Lithosphere with rift"); + { + LAB_isotherm = prm.get_double ("LAB isotherm temperature"); + T0 = prm.get_double ("Surface temperature"); + use_compensation_depth = prm.get_bool ("Use temperature compensation depth"); + compensation_depth = prm.get_double ("Temperature compensation depth"); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + + + // Retrieve the indices of the fields that represent the lithospheric layers. + AssertThrow(this->introspection().compositional_name_exists("upper"),ExcMessage("We need a compositional field called 'upper' representing the upper crust.")); + AssertThrow(this->introspection().compositional_name_exists("lower"),ExcMessage("We need a compositional field called 'lower' representing the lower crust.")); + AssertThrow(this->introspection().compositional_name_exists("mantle_L"),ExcMessage("We need a compositional field called 'mantle_L' representing the lithospheric part of the mantle.")); + + // For now, we assume a 3-layer system with an upper crust, lower crust and lithospheric mantle + const unsigned int id_upper = this->introspection().compositional_index_for_name("upper"); + const unsigned int id_lower = this->introspection().compositional_index_for_name("lower"); + const unsigned int id_mantle_L = this->introspection().compositional_index_for_name("mantle_L"); + + // Retrieve other material properties set in different sections such that there + // is no need to set them twice. + + prm.enter_subsection("Heating model"); + { + prm.enter_subsection("Compositional heating"); + { + // The heating model compositional heating prefixes an entry for the background material + const std::vector temp_heat_productivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Compositional heating values"))), + n_fields+1, + "Compositional heating values"); + // This sets the heat productivity in W/m3 units + heat_productivities.push_back(temp_heat_productivities[id_upper+1]); + heat_productivities.push_back(temp_heat_productivities[id_lower+1]); + heat_productivities.push_back(temp_heat_productivities[id_mantle_L+1]); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + + prm.enter_subsection ("Compositional fields"); + { + list_of_composition_names = Utilities::split_string_list (prm.get("Names of fields")); + } + prm.leave_subsection(); + + prm.enter_subsection("Material model"); + { + prm.enter_subsection("Visco Plastic"); + { + const std::vector temp_thermal_diffusivities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Thermal diffusivities"))), + n_fields+1, + "Thermal diffusivities"); + const std::vector temp_heat_capacities = Utilities::possibly_extend_from_1_to_N (Utilities::string_to_double(Utilities::split_string_list(prm.get("Heat capacities"))), + n_fields+1, + "Heat capacities"); + + n_phases_for_each_composition = std::make_unique>(); + const std::vector temp_densities = Utilities::parse_map_to_double_array (prm.get("Densities"), + list_of_composition_names, + /*has_background_field=*/true, + "Densities", + true, + n_phases_for_each_composition); + + // Assemble a list of phase densities for each composition. + // Add 1 for background material. + std::vector> densities_per_composition(this->n_compositional_fields()+1); + unsigned int counter = 0; + for (unsigned int i = 0; i < (*n_phases_for_each_composition).size(); ++i) + { + for (unsigned int j = 0; j < (*n_phases_for_each_composition)[i]; ++j) + { + densities_per_composition[i].push_back(temp_densities[counter]); + ++counter; + } + } + + densities.push_back(densities_per_composition[id_upper+1][0]); + densities.push_back(densities_per_composition[id_lower+1][0]); + densities.push_back(densities_per_composition[id_mantle_L+1][0]); + + // Thermal diffusivity kappa = k/(rho*cp), so thermal conducitivity k = kappa*rho*cp + conductivities.push_back(temp_thermal_diffusivities[id_upper+1] * densities[0] * temp_heat_capacities[id_upper+1]); + conductivities.push_back(temp_thermal_diffusivities[id_lower+1] * densities[1] * temp_heat_capacities[id_lower+1]); + conductivities.push_back(temp_thermal_diffusivities[id_mantle_L+1] * densities[2] * temp_heat_capacities[id_mantle_L+1]); + + // To obtain the radioactive heating rate in W/kg, we divide the volumetric heating rate by density + AssertThrow(heat_productivities.size() == 3 && densities.size() == 3 && conductivities.size() == 3, + ExcMessage("The entries for density, conductivity and heat production do not match with the expected number of layers (3).")); + + for (unsigned int i = 0; i<3; ++i) + heat_productivities[i] /= densities[i]; + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + } + + +} + + + + +// explicit instantiations +namespace aspect +{ + namespace InitialTemperature + { + ASPECT_REGISTER_INITIAL_TEMPERATURE_MODEL(LithosphereRift, + "lithosphere with rift", + "This is a temperature initial condition that " + "computes a continental geotherm based on the " + "conductive equations of Turcotte and Schubert Ch. 4.6. " + "This plugin only works with the corresponding composition " + "initial conditions plugin 'lithosphere with rift'. " + "The geotherm can be computed for any number of crustal " + "layers, for each of which a density, heat production and thermal " + "conductivity should be supplied. " + "Make sure the top and bottom temperatures of the lithosphere " + "agree with temperatures set in for example the temperature " + "boundary conditions. " + "The thickness of the lithospheric layers is adapted for distance " + "to the rift specified with the InitialComposition plugin `lithosphere " + "with rift'. Additional variations in lithospheric thickness can be " + "added through specifying polygonal areas of different thicknesses. ") + } +} diff --git a/cookbooks/initial_conditions_continental_rift/lithosphere_rift_topo.cc b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_topo.cc new file mode 100644 index 00000000000..b13da041415 --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_topo.cc @@ -0,0 +1,256 @@ +/* + Copyright (C) 2016 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 + . +*/ + + +#include "lithosphere_rift_topo.h" +#include "lithosphere_rift_IC.h" +#include +#include +#include +#include + +namespace aspect +{ + namespace InitialTopographyModel + { + template + void + LithosphereRift:: + initialize () + { + // Compute the maximum topography amplitude based on isostasy. + // Assume the reference density is representative for each layer (despite temperature dependence) + + // For now, we assume a 3-layer system with an upper crust, lower crust and lithospheric mantle + const unsigned int id_upper = this->introspection().compositional_index_for_name("upper"); + const unsigned int id_lower = this->introspection().compositional_index_for_name("lower"); + const unsigned int id_mantle_L = this->introspection().compositional_index_for_name("mantle_L"); + + // Assemble a list of phase densities for each composition. + // Add 1 for background material. + std::vector> densities_per_composition(this->n_compositional_fields()+1); + unsigned int counter = 0; + for (unsigned int i = 0; i < (*n_phases_for_each_composition).size(); ++i) + { + for (unsigned int j = 0; j < (*n_phases_for_each_composition)[i]; ++j) + { + densities_per_composition[i].push_back(temp_densities[counter]); + ++counter; + } + } + + // Get the relevant densities for the lithosphere. + // We take the reference density of the first phase. + densities.push_back(densities_per_composition[0][0]); + densities.push_back(densities_per_composition[id_upper+1][0]); + densities.push_back(densities_per_composition[id_lower+1][0]); + densities.push_back(densities_per_composition[id_mantle_L+1][0]); + + // The reference column + ref_rgh = 0; + // Assume constant gravity magnitude, so ignore + for (unsigned int l=0; l<3; ++l) + ref_rgh += densities[l+1] * reference_thicknesses[l]; + + // The total lithosphere thickness + const double sum_thicknesses = std::accumulate(reference_thicknesses.begin(), reference_thicknesses.end(), 0); + + // The column at the rift center + double rift_rgh = 0; + rift_thicknesses = reference_thicknesses; + for (unsigned int l=0; l polygon_rgh(n_polygons); + std::vector sum_polygon_thicknesses(n_polygons); + double max_sum_polygon_thicknesses = 0.; + for (unsigned int i_polygons=0; i_polygonsget_pcout() << " Maximum initial topography of rift: " << topo_rift_amplitude << " m" << std::endl; + this->get_pcout() << " Maximum initial topography of polygon: " << topo_polygon_amplitude << " m" << std::endl; + } + + + template + double + LithosphereRift:: + value (const Point &position) const + { + // The simulator only keeps the initial conditions around for + // the first time step. As a consequence, we have to save a + // shared pointer to that object ourselves the first time we get + // here. + if (initial_composition_manager == nullptr) + const_cast>&>(initial_composition_manager) = this->get_initial_composition_manager_pointer(); + + // Check that the required initial composition model is used + // We have to do it here instead of in initialize() because + // the names are not available upon initialization of the + // initial topography model yet. + const std::vector active_initial_composition_models = initial_composition_manager->get_active_initial_composition_names(); + AssertThrow(initial_composition_manager->template has_matching_initial_composition_model>(), + ExcMessage("The 'lithosphere with rift' initial topography plugin requires the 'lithosphere with rift' initial composition plugin.")); + + // When cartesian, position contains x(,y); when spherical, position contains lon(,lat) (in degrees). + // Turn into a Point + Point surface_position; + for (unsigned int d=0; d &ic = initial_composition_manager->template get_matching_initial_composition_model>(); + + // Compute the topography based on distance to the rift and distance to the polygon + std::vector local_thicknesses = ic.compute_local_thicknesses(surface_position); + + // The local lithospheric column + double local_rgh = 0; + for (unsigned int l=0; l<3; ++l) + local_rgh += densities[l+1] * local_thicknesses[l]; + // The total local lithosphere thickness + const double sum_local_thicknesses = std::accumulate(local_thicknesses.begin(), local_thicknesses.end(),0); + local_rgh += (compensation_depth - sum_local_thicknesses) * densities[0]; + + return (ref_rgh - local_rgh) / densities[0]; + } + + + template + double + LithosphereRift:: + max_topography () const + { + return maximum_topography; + } + + + template + void + LithosphereRift:: + declare_parameters (ParameterHandler &) + { + } + + + + template + void + LithosphereRift::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection ("Initial composition model"); + { + prm.enter_subsection("Lithosphere with rift"); + { + A_rift = Utilities::possibly_extend_from_1_to_N(Utilities::string_to_double(Utilities::split_string_list(prm.get("Amplitude of Gaussian rift geometry"))), + 3, + "Amplitude of Gaussian rift geometry"); + reference_thicknesses = Utilities::possibly_extend_from_1_to_N(Utilities::string_to_double(Utilities::split_string_list(prm.get("Layer thicknesses"))), + 3, + "Layer thicknesses"); + // Split the string into the separate polygons + const std::vector temp_thicknesses = Utilities::split_string_list(prm.get("Lithospheric polygon layer thicknesses"),';'); + const unsigned int n_polygons = temp_thicknesses.size(); + polygon_thicknesses.resize(n_polygons); + for (unsigned int i_polygons = 0; i_polygons < n_polygons; ++i_polygons) + { + polygon_thicknesses[i_polygons] = Utilities::string_to_double(Utilities::split_string_list(temp_thicknesses[i_polygons],',')); + AssertThrow(polygon_thicknesses[i_polygons].size()==3, ExcMessage ("The number of layer thicknesses should be equal to 3.")); + } + } + prm.leave_subsection(); + } + prm.leave_subsection(); + + prm.enter_subsection ("Compositional fields"); + { + list_of_composition_names = Utilities::split_string_list (prm.get("Names of fields")); + } + prm.leave_subsection(); + + prm.enter_subsection("Material model"); + { + prm.enter_subsection("Visco Plastic"); + { + n_phases_for_each_composition = std::make_unique>(); + temp_densities = Utilities::parse_map_to_double_array (prm.get("Densities"), + list_of_composition_names, + /*has_background_field=*/true, + "Densities", + true, + n_phases_for_each_composition); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + } +} + +// explicit instantiations +namespace aspect +{ + namespace InitialTopographyModel + { + ASPECT_REGISTER_INITIAL_TOPOGRAPHY_MODEL(LithosphereRift, + "lithosphere with rift", + "An initial topography model that defines the initial topography " + "based on isostasy. It takes into account lithospheric thickness " + "variations as specified in the InitialComposition model lithosphere " + "with rift' as should only be used in conjunction with this model. ") + } +} diff --git a/cookbooks/initial_conditions_continental_rift/lithosphere_rift_topo.h b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_topo.h new file mode 100644 index 00000000000..7020543fd8f --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/lithosphere_rift_topo.h @@ -0,0 +1,160 @@ +/* + Copyright (C) 2024 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 + . +*/ + + +#ifndef _aspect_geometry_model__initial_topography_lithosphere_rift_h +#define _aspect_geometry_model__initial_topography_lithosphere_rift_h + +#include +#include + +namespace aspect +{ + namespace InitialTopographyModel + { + using namespace dealii; + + /** + * A class that describes an initial topography for the geometry model, + * based on isostatically balancing the thinning of lithospheric layers + * around a rift axis. + */ + template + class LithosphereRift : public Interface, + public SimulatorAccess + { + public: + /** + * Initialization function. This function is called once at the + * beginning of the program. Checks preconditions. + */ + void + initialize () override; + + /** + * Return the value of the topography for a point. + */ + virtual + double value (const Point &p) const override; + + /** + * Return the maximum value of the elevation. + */ + virtual + double max_topography () const override; + + /** + * Declare the parameters this class takes through input files. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + */ + virtual + void + parse_parameters (ParameterHandler &prm) override; + + private: + + /** + * The maximum amplitude of the Gaussian distribution of the thinning/thickening + * of the lithospheric thicknesses with distance from the rift axis. + * The amplitude should have values between -1 and 1, where positive + * numbers represent a reduction in thickness and negative numbers an increase. + * For example, values of 0.25, 0, 0 reduce the reference thickness of the + * upper crust by 25%, while the lower crust and mantle lithosphere are + * untouched. + */ + std::vector A_rift; + + /** + * The maximum amplitude of the Gaussian distribution of the topography around the rift. + */ + double topo_rift_amplitude; + + /** + * The product of density, gravitational acceleration constant and thickness of the + * reference lithospheric column used in computing the topography based on isostasy. + */ + double ref_rgh; + + /** + * The thickness/depth of the reference lithospheric column used in computing the topography + * based on isostasy. + */ + double compensation_depth; + + /** + * The maximum amplitude of the topography of the polygon area. + */ + double topo_polygon_amplitude; + + /** + * Vector for field densities. + */ + std::vector densities; + std::vector temp_densities; + + /** + * Vector for the reference field thicknesses away from the rift. + */ + std::vector reference_thicknesses; + + /** + * Vector for the rift thicknesses at the center. + */ + std::vector rift_thicknesses; + + /** + * Vector for the polygon thicknesses at the center. + */ + std::vector> polygon_thicknesses; + + /** + * The maximum topography in this model. + */ + double maximum_topography; + + /** + * A shared pointer to the initial composition object + * that ensures that the current object can continue + * to access the initial composition object beyond the + * first time step. + */ + std::shared_ptr> initial_composition_manager; + + /** + * Vector containing the number of phases for each composition. + */ + std::unique_ptr> n_phases_for_each_composition; + + /** + * Vector containing the names of the compositional fields. + */ + std::vector list_of_composition_names; + }; + } +} + + +#endif diff --git a/cookbooks/initial_conditions_continental_rift/rift_box_initial_plastic_strain.cc b/cookbooks/initial_conditions_continental_rift/rift_box_initial_plastic_strain.cc new file mode 100644 index 00000000000..59624a6b47f --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/rift_box_initial_plastic_strain.cc @@ -0,0 +1,339 @@ +/* + Copyright (C) 2011 - 2016 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 doc/COPYING. If not see + . + */ + + +#include "rift_box_initial_plastic_strain.h" +#include +#include +#include +#include + +namespace aspect +{ + namespace InitialComposition + { + template + RiftBoxInitialPlasticStrain::RiftBoxInitialPlasticStrain () + {} + + template + void + RiftBoxInitialPlasticStrain::initialize () + { + AssertThrow(Plugins::plugin_type_matches>(this->get_material_model()), + ExcMessage("This initial condition only makes sense in combination with the visco_plastic material model.")); + + + // From shear_bands.cc + Point extents_max; + TableIndices size_idx; + for (unsigned int d=0; d white_noise; + white_noise.TableBase::reinit(size_idx); + + // Max of each direction (m) + extents_max = extents + origin; + + // Add initial topography to vertical direction + extents_max[dim-1] += maximum_initial_topography; + + for (unsigned int d=0; d idx; + + for (unsigned int i=0; i (grid_extents, + grid_intervals, + white_noise); + } + + template + double + RiftBoxInitialPlasticStrain:: + initial_composition (const Point &position, const unsigned int n_comp) const + { + // If n_comp does not represent the strain field, + // return 0 right away. + if (n_comp != strain_composition_number) + return 0.0; + + // Initiate distance with large value + double distance_to_rift_axis = 1e23; + double temp_distance = 0; + + const Point natural_coords = position; + + // Loop over all line segments + for (unsigned int i_segments = 0; i_segments < point_list.size(); ++i_segments) + { + if (dim == 2) + temp_distance = std::abs(natural_coords[0]-point_list[i_segments][0][0]); + else + { + // Get the surface coordinates by dropping the last coordinate + const Point<2> surface_position = Point<2>(natural_coords[0],natural_coords[1]); + temp_distance = std::abs(Utilities::distance_to_line(point_list[i_segments], surface_position)); + } + + // Get the minimum distance + distance_to_rift_axis = std::min(distance_to_rift_axis, temp_distance); + } + + // Smoothing of noise with height (depth with respect to the unperturbed surface) + const double depth_smoothing = 0.5 * (1.0 + std::tanh((position[dim-1] - strain_height) / strain_halfwidth)); + // Smoothing of noise with lateral distance to the rift axis + const double noise_amplitude = A * std::exp((-std::pow(distance_to_rift_axis,2)/(2.0*std::pow(sigma,2)))) * depth_smoothing; + + // Add randomness + return noise_amplitude * interpolate_noise->value(natural_coords); + } + + template + void + RiftBoxInitialPlasticStrain::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Initial composition model"); + { + prm.enter_subsection("Rift box initial plastic strain"); + { + prm.declare_entry ("Random number generator seed", "0", + Patterns::Integer (0), + "The value of the seed used for the random number generator. " + "Units: none."); + prm.declare_entry ("Standard deviation of Gaussian noise amplitude distribution", "20000", + Patterns::Double (0), + "The standard deviation of the Gaussian distribution of the amplitude of the strain noise. " + "Note that this parameter is taken to be the same for all rift segments. " + "Units: $m$ or degrees."); + prm.declare_entry ("Maximum amplitude of Gaussian noise amplitude distribution", "0.2", + Patterns::Double (0), + "The amplitude of the Gaussian distribution of the amplitude of the strain noise. " + "Note that this parameter is taken to be the same for all rift segments. " + "Units: none."); + prm.declare_entry ("Depth around which Gaussian noise is smoothed out", "40000", + Patterns::Double (0), + "The depth around which smoothing out of the strain noise starts with a hyperbolic tangent. " + "Note that this parameter is taken to be the same for all rift segments. Also, depth " + "is with respect to the unperturbed surface of the box domain. " + "Units: $m$."); + prm.declare_entry ("Halfwidth with which Gaussian noise is smoothed out in depth", "40000", + Patterns::Double (0), + "The halfwidth with which smoothing out of the strain noise is done with a hyperbolic tangent. " + "Note that this parameter is taken to be the same for all rift segments. " + "Units: $m$."); + prm.declare_entry ("Rift axis line segments", + "", + Patterns::Anything(), + "Set the line segments that represent the rift axis. In 3d each segment is made up of " + "two points that represent horizontal coordinates (x,y). " + "The exact format for the point list describing the segments is " + "\"x1,y1>x2,y2;x2,y2>x3,y3;x4,y4>x5,y5\". In 2d, a segment is made up by 1 horizontal " + "x coordinate: \"x1;x2;x3\". Note that the segments can be connected " + "or isolated. Units: $m$."); + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + + template + void + RiftBoxInitialPlasticStrain::parse_parameters (ParameterHandler &prm) + { + // Check that there is a compositional field called strain and retrieve its index + AssertThrow(this->introspection().compositional_name_exists("plastic_strain"), + ExcMessage("This plugin requires a compositional field named plastic_strain. ")); + strain_composition_number = this->introspection().compositional_index_for_name("plastic_strain"); + + AssertThrow(Plugins::plugin_type_matches>(this->get_geometry_model()), + ExcMessage("This initial condition only makes sense in combination with the box geometry model.")); + + // Get the number of mesh cells in each direction + prm.enter_subsection("Geometry model"); + { + prm.enter_subsection("Box"); + { + origin[0] = prm.get_double ("Box origin X coordinate"); + extents[0] = prm.get_double ("X extent"); + repetitions[0] = prm.get_integer ("X repetitions"); + origin[1] = prm.get_double ("Box origin Y coordinate"); + extents[1] = prm.get_double ("Y extent"); + repetitions[1] = prm.get_integer ("Y repetitions"); + if (dim >= 3) + { + origin[dim-1] = prm.get_double ("Box origin Z coordinate"); + extents[dim-1] = prm.get_double ("Z extent"); + repetitions[dim-1] = prm.get_integer ("Z repetitions"); + } + } + prm.leave_subsection(); + } + prm.leave_subsection(); + + unsigned int refinement; + prm.enter_subsection("Mesh refinement"); + { + refinement = prm.get_integer("Initial adaptive refinement") + prm.get_integer("Initial global refinement"); + } + prm.leave_subsection(); + + // The grid intervals equal the number of mesh cells. + // refinement, repetitions and grid_intervals are integers, + // but pow returns a double, so convert to int. + for (unsigned int d = 0; d < dim; ++d) + grid_intervals[d] = repetitions[d] * int(std::pow(2, refinement)); + + // Both initial mesh deformation and initial topography models + // can distort the initial mesh such that the vertical coordinate + // of some locations is larger than the initial, unperturbed, + // vertical extent of the model domain. Since not all initial topography + // models (correctly) return a maximum topography, and + // initial mesh deformation models do not at all, we here assume + // a maximal initial topography of 10 km. This is added to the + // vertical extent. + // + // In the vertical direction, we add enough grid cells + // to cover 10 km at the same resolution. Round up the number of cells. + grid_intervals[dim-1] += static_cast(std::ceil(maximum_initial_topography / (extents[dim-1] / grid_intervals[dim-1]))); + + prm.enter_subsection("Initial composition model"); + { + prm.enter_subsection("Rift box initial plastic strain"); + sigma = prm.get_double ("Standard deviation of Gaussian noise amplitude distribution"); + A = prm.get_double ("Maximum amplitude of Gaussian noise amplitude distribution"); + seed = prm.get_integer ("Random number generator seed"); + strain_height = origin[dim-1] + extents[dim-1] - prm.get_double ("Depth around which Gaussian noise is smoothed out"); + strain_halfwidth = prm.get_double ("Halfwidth with which Gaussian noise is smoothed out in depth"); + + // Read in the string of rift segments + const std::string temp_all_segments = prm.get("Rift axis line segments"); + // Split the string into segment strings + const std::vector temp_segments = Utilities::split_string_list(temp_all_segments,';'); + // The number of segments, each consisting of a begin and an end point in 3d and one point in 3d + const unsigned int n_temp_segments = temp_segments.size(); + point_list.resize(n_temp_segments); + + // Loop over the segments to extract the points + for (unsigned int i_segment = 0; i_segment < n_temp_segments; i_segment++) + { + // In 3d a line segment consists of 2 points, + // in 2d of only 1 (ridge axis orthogonal to x and y). + // Also, a 3d point has 2 coordinates (x and y), + // a 2d point only 1 (x). + const std::vector temp_segment = Utilities::split_string_list(temp_segments[i_segment],'>'); + + if (dim == 3) + { + // const std::vector temp_segment = Utilities::split_string_list(temp_segments[i_segment],'>'); + AssertThrow(temp_segment.size() == 2,ExcMessage ("The given coordinate '" + temp_segment[i_segment] + "' is not correct. " + "It should only contain 2 parts: " + "the two points of the segment, separated by a '>'.")); + + } + else + { + // Add the point to the list of points for this segment + // As we're in 2d all segments correspond to 1 point consisting of 1 coordinate + AssertThrow(temp_segment.size() == 1,ExcMessage ("The given coordinate '" + temp_segment[i_segment] + "' is not correct. " + "It should only contain 1 part: " + "the point representing the rift axis.")); + } + + // Loop over the 2 points of each segment + for (unsigned int i_points = 0; i_points < dim-1; i_points++) + { + std::vector temp_point = Utilities::string_to_double(Utilities::split_string_list(temp_segment[i_points],',')); + if (dim == 3) + { + AssertThrow(temp_point.size() == 2,ExcMessage ("The given coordinates of segment '" + temp_segment[i_points] + "' are not correct. " + "It should only contain 2 parts: " + "the x and y coordinates of the segment begin/end point, separated by a ','.")); + } + else + { + AssertThrow(temp_point.size() == 1,ExcMessage ("The given coordinates of segment '" + temp_segment[i_points] + "' are not correct. " + "It should only contain 1 part: " + "the one coordinate of the segment end point.")); + } + + // Add the point to the list of points for this segment + point_list[i_segment][i_points][0] = temp_point[0]; + point_list[i_segment][i_points][1] = temp_point[dim-2]; + } + } + prm.leave_subsection(); + } + prm.leave_subsection(); + } + + } +} + +// explicit instantiations +namespace aspect +{ + namespace InitialComposition + { + ASPECT_REGISTER_INITIAL_COMPOSITION_MODEL(RiftBoxInitialPlasticStrain, + "rift box initial plastic strain", + "Specify the plastic strain initial compositional field value based on " + "the distance to a list of line segments " + "and the user-defined Gaussian distribution around these segments, " + "combined with random noise. " + "A maximum initial topography of 10 km can be accommodated. Initial topography that " + "is higher than 10 km will lead to plastic strain values set to the strain value at 10 km." + "This plugin only works in combination with the box geometry model and the visco_plastic " + "material model.") + } +} diff --git a/cookbooks/initial_conditions_continental_rift/rift_box_initial_plastic_strain.h b/cookbooks/initial_conditions_continental_rift/rift_box_initial_plastic_strain.h new file mode 100644 index 00000000000..831dc9ab86e --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/rift_box_initial_plastic_strain.h @@ -0,0 +1,162 @@ +/* + Copyright (C) 2011 - 2016 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 doc/COPYING. If not see + . +*/ + + +#ifndef _aspect_initial_composition_rift_box_initial_plastic_strain_h +#define _aspect_initial_composition_rift_box_initial_plastic_strain_h + +#include +#include + +#include + +namespace aspect +{ + namespace InitialComposition + { + using namespace dealii; + + /** + * A class that implements initial conditions for the compositional fields + * based on a functional description provided in the input file. + * + * @ingroup InitialCompositionModels + */ + template + class RiftBoxInitialPlasticStrain : public Interface, + public SimulatorAccess + { + public: + /** + * Constructor. + */ + RiftBoxInitialPlasticStrain (); + + /** + * Initialization function. + */ + void + initialize () override; + + /** + * Return the initial composition as a function of position and number + * of compositional field. The composition varies as a Gaussian distribution + * around a user-defined set of line-segments. + */ + virtual + double initial_composition (const Point &position, const unsigned int n_comp) const override; + + /** + * Declare the parameters this class takes through input files. The + * default implementation of this function does not describe any + * parameters. Consequently, derived classes do not have to overload + * this function if they do not take any runtime parameters. + */ + static + void + declare_parameters (ParameterHandler &prm); + + /** + * Read the parameters this class declares from the parameter file. + * The default implementation of this function does not read any + * parameters. Consequently, derived classes do not have to overload + * this function if they do not take any runtime parameters. + */ + virtual + void + parse_parameters(ParameterHandler &prm) override; + + private: + + /** + * The number of the compositional field representing the strain. + */ + unsigned int strain_composition_number; + + /** + * The value of the seed for the random number generator + */ + unsigned int seed; + + /** + * The maximum amplitude of the Gaussian amplitude of the noise. + */ + double A; + + /** + * The standard deviation of the Gaussian amplitude of the noise. + */ + double sigma; + + /** + * The height around and halfwidth with which the noise is smoothed out. + */ + double strain_height; + double strain_halfwidth; + + /** + * The list of line segments consisting of two 2d coordinates per segment (even in 2d). + * The segments represent the rift axis. + */ + std::vector,2 >> point_list; + + /** + * A table with random noise for the + * second invariant of the strain. + */ + std::array grid_intervals; + Functions::InterpolatedUniformGridData *interpolate_noise; + + /** + * A pointer to the initial topography model. + */ + InitialTopographyModel::Interface *topo_model; + + /** + * The min and max of the strain grid in each direction + */ + std::array,dim> grid_extents; + + /** + * Origin of the box in x, y, and z (in 3d) coordinates. + */ + Point origin; + + /** + * Extent of the domain in x-, y-, and z-direction (in 3d). + */ + Point extents; + + /** + * The number of cells in each coordinate direction (x, y, z). + */ + std::array repetitions; + + /** + * The height we add to vertical domain extent + * to deal with any initial topography. + */ + const double maximum_initial_topography = 10000.; + }; + } +} + + +#endif diff --git a/cookbooks/initial_conditions_continental_rift/test.prm b/cookbooks/initial_conditions_continental_rift/test.prm new file mode 100644 index 00000000000..ac93ccba0cd --- /dev/null +++ b/cookbooks/initial_conditions_continental_rift/test.prm @@ -0,0 +1,251 @@ +#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 = 3.5e3 +set Use years in output instead of seconds = true +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-6 +set Max nonlinear iterations in pre-refinement = 0 + + +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 + subsection Initial topography model + set Model name = lithosphere with rift + end +end + +subsection Initial temperature model + set List of model names = adiabatic + subsection Adiabatic + 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 = function + + 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 Function expression = 0; \ + 0; \ + 0; \ + 0; \ + 0; \ + 0; \ + 0; \ + 0; \ + 0; \ + 0 + end +end + +subsection Boundary temperature model + set List of model names = box + set Fixed temperature boundary indicators = bottom, top + subsection Box + set Top temperature = 293 + set Bottom temperature = 1700 # Unrealistically high, so that it is always taken from initial temperature plugin + end +end + + +subsection Mesh deformation + set Mesh deformation boundary indicators = top: free surface +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 Tangential velocity boundary indicators = left, right, bottom +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 = 2 + set Initial adaptive refinement = 0 #3 + 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 + +end diff --git a/doc/sphinx/user/cookbooks/geophysical-setups.md b/doc/sphinx/user/cookbooks/geophysical-setups.md index 6146b3d9cfc..bf113543c55 100644 --- a/doc/sphinx/user/cookbooks/geophysical-setups.md +++ b/doc/sphinx/user/cookbooks/geophysical-setups.md @@ -102,6 +102,7 @@ cookbooks/multicomponent_steinberger/doc/multicomponent_steinberger.md cookbooks/morency_doin_2004/doc/morency_doin_2004.md cookbooks/crustal_deformation/doc/crustal_deformation.md cookbooks/continental_extension/doc/continental_extension.md +cookbooks/initial_conditions_continental_rift/doc/initial_conditions_continental_rift.md cookbooks/inner_core_convection/doc/inner_core_convection.md cookbooks/lower_crustal_flow/doc/lower_crustal_flow.md cookbooks/global_melt/doc/global_melt.md