Skip to content

Commit

Permalink
Improve continent cookbook
Browse files Browse the repository at this point in the history
  • Loading branch information
@gassmoeller @jdannberg
gassmoeller authored and jdannberg committed Jun 12, 2024
1 parent 805fdc3 commit 9814c9f
Show file tree
Hide file tree
Showing 2 changed files with 47 additions and 38 deletions.
28 changes: 14 additions & 14 deletions ...antle_convection_with_continents_in_annulus/doc/mantle_convection_in_annulus.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,33 +4,33 @@
*This section was contributed by Cedric Thieulot and Erik van der Wiel.*

The setup for this experiment originates in {cite:t}`vanderwiel:etal:2024a`.
The domain is an annulus with an inner radius of 3480km and an outer radius of 6370km.
The Isentropic Compression Approximation (ICA) is used {cite:t}`gassmoller:etal:2020`, which is
The domain is an annulus with an inner radius of 3480 km and an outer radius of 6370 km.
The Isentropic Compression Approximation (ICA) is used {cite}`gassmoller:etal:2020`, which is
the default approximation for compressible flow in ASPECT.

Flow in the model is governed by the visco-plastic flow equations
that describe dislocation and diffusion creep and we use the Drucker-Prager
yield criterion to limit viscous stresses {cite:t}`glerum:etal:2018`.
The grain-size in the diffusion creep is assumed constant and therefore
yield criterion to limit viscous stresses {cite}`glerum:etal:2018`.
The grain size in the diffusion creep is assumed constant and therefore
incorporated in the pre-factor A.
Note that in the lower mantle flow is based on diffusion creep only.

Both the inner and outer boundaries are free-slip boundaries so
there is no external kinematic forcing on the model. The resulting
existing rotational null space is removed by imposing no-net-rotation of
the mantle. The boundaries have fixed temperatures of 3700K and 300K.
the mantle. The boundaries have fixed temperatures of 3700 K and 300 K.
Within the domain two compositional fields are defined: mantle and
continent. The mantle domain consists of three different regions separated
by the two major phase changes occurring at ~410- and ~660 km
depths. We use the geodynamic World Builder {cite:t}`Fraters2019c` to
depths. We use the geodynamic World Builder {cite}`Fraters2019c` to
set the initial temperature and compositional field distribution
(see {numref}`fig:init_visc`, {numref}`fig:init_compositions` and {numref}`fig:init_temperature`).


```{figure-md} fig:init_visc
<img src="viscosity.*" width="90%" />
the initial viscosity field.
Initial viscosity field.
```

```{figure-md} fig:init_compositions
Expand All @@ -42,23 +42,23 @@ The initial position of the three continents.
```{figure-md} fig:init_temperature
<img src="temperature.*" width="90%" />
The initial temperature field.
Initial temperature field.
```

Similar to {cite:t}`ulvrova:etal:2019`, we use viscous rafts as `continents' to aid with
Similar to {cite:t}`ulvrova:etal:2019`, we use viscous rafts as "continents" to aid with
modelling one-sided subduction systems. We use three continents of
5000, 5000 and 3000 km length, covering roughly 30% of the model
surface. To avoid subduction of the continents, we assign the viscous
rafts a reference density of 2916 kg/m3, which corresponds to a 400
kg/m3 density contrast with the reference mantle density.
rafts a reference density of 2916 kg/m<sup>3</sup>, which corresponds to a 400
kg/m<sup>3</sup> density contrast with the reference mantle density.
Furthermore, to avoid deformation of the continents, we assign
them a viscosity that equals the maximum cut-off viscosity in the model,
i.e. 6e24 Pa s. We also configure three subduction zones
in the initial set-up to seed the model with initial negative
buoyancy.

As the goal of the publication is to investigate
how average slab sinking rates relate to the vigour of mantle convection and mixing,
As the goal of the cookbook (and the corresponding publication in {cite:t}`vanderwiel:etal:2024a`) is to investigate
how average slab sinking rates relate to the vigor of mantle convection and mixing,
various rheologies are considered in the lower mantle, as
shown in {numref}`fig:visc-profiles`.
We will focus on the 'R' (reference) profile in this cookbook (green line).
Expand All @@ -83,7 +83,7 @@ Results of the reference model (R) after 500 Ma of mantle convection simulation.
temperature (b) within the modelled domain, showing five slabs actively subducting below the continents (pink).
(c) Solid lines indicate average surface velocity (blue)
and mobility M (red) of model R as well as their dashed time-averages Vsurf=2.1 cm/a and M = 2.218.
Mobility is defined as the ratio of rms surface velocity to rms velocity averaged over the entire 3D domain{cite:t}`tackley:2000`.
Mobility is defined as the ratio of rms surface velocity to rms velocity averaged over the entire 3D domain{cite}`tackley:2000`.
(d) Radial velocity of all tracers defined as
slabs plotted at their position in the model where blue indicates sinking slabs and red slowly rising. Shown tracers are 300 K colder than the radial averaged
temperature at similar depth and automatically obtained (see methods section of the publication).
Expand Down
57 changes: 33 additions & 24 deletions cookbooks/mantle_convection_with_continents_in_annulus/modelR.prm
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# input file for model R of van der wiel et al, 2024.
# input file for model R of van der Wiel et al., 2024.

set Dimension = 2
set Use years in output instead of seconds = true
Expand All @@ -12,36 +12,41 @@ set Nonlinear solver tolerance = 1e-4
set Timing output frequency = 10

subsection Solver parameters
subsection Stokes solver parameters
set Number of cheap Stokes solver steps = 0
set Linear solver tolerance = 1e-6
set Maximum number of expensive Stokes solver steps = 5000
end
subsection Stokes solver parameters
set Stokes solver type = block GMG
set Number of cheap Stokes solver steps = 5000
set Linear solver tolerance = 1e-6
end
end

subsection Formulation
set Formulation = isentropic compression
end

subsection Compositional fields
set Number of fields = 2
set Names of fields = 1_mantle, 3_continent
set Number of fields = 2
set Names of fields = 1_mantle, 3_continent
set Compositional field methods = particles, particles
set Mapped particle properties = 1_mantle: initial 1_mantle, 3_continent: initial 3_continent
end

subsection Nullspace removal
set Remove nullspace = net rotation
end

subsection Material model

set Model name = compositing
set Material averaging = harmonic average only viscosity

subsection Visco Plastic
set Reference temperature = 1600
set Minimum viscosity = 1e20
set Maximum viscosity = 6e24
set Viscous flow law = composite
set Minimum strain rate = 1e-22
set Minimum viscosity = 1e20
set Maximum viscosity = 6e24
set Viscous flow law = composite
set Minimum strain rate = 1e-22
set Reference strain rate = 1e-15
set Heat capacities = background: 1250, 1_mantle: 1250| 1250 | 1250 , 3_continent:1250
set Heat capacities = background: 1250, 1_mantle: 1250 | 1250 | 1250 , 3_continent: 1250

set Prefactors for dislocation creep = background:6.51e-16 , 1_mantle: 6.51e-16| 8.51e-16 | 6.51e-16 , 3_continent:6.51e-28
set Stress exponents for dislocation creep = background:1 , 1_mantle: 3 | 3 | 1 , 3_continent: 1
Expand All @@ -63,7 +68,7 @@ subsection Material model
set Phase transition widths = 1_mantle: 50000 | 50000

set Use adiabatic pressure in creep viscosity = true
set Constant viscosity prefactors = 1, 1, 1
set Constant viscosity prefactors = 1, 1, 1
end

subsection Latent heat
Expand All @@ -81,8 +86,8 @@ subsection Material model
# for T=T1.
set Phase transition depths = 1_mantle: 410000 | 660000
set Phase transition temperatures = 1_mantle: 1700 | 1900
set Phase transition Clapeyron slopes = 1_mantle: 3e6 | -2.5e6
set Phase transition widths = 1_mantle: 50000 | 50000
set Phase transition Clapeyron slopes = 1_mantle: 3e6 | -2.5e6
set Phase transition widths = 1_mantle: 50000 | 50000
set Reference specific heat = 1250
set Thermal conductivity = 6
set Thermal expansion coefficient = 3e-5
Expand Down Expand Up @@ -121,7 +126,7 @@ subsection Boundary velocity model
end

subsection Heating model
set List of model names = latent heat, adiabatic heating, constant heating, shear heating
set List of model names = latent heat, adiabatic heating, constant heating, shear heating
subsection Adiabatic heating
set Use simplified adiabatic heating = true
end
Expand Down Expand Up @@ -178,26 +183,30 @@ subsection Postprocess
set Number of grouped files = 1
set Output format = vtu
set Time between graphical output = 1e6
set List of output variables = density, viscosity, strain rate, stress, temperature anomaly, spherical velocity components,adiabat, heat flux map, heating, vertical heat flux
set List of output variables = density, viscosity, strain rate, stress, temperature anomaly, spherical velocity components, adiabat, heat flux map, heating, vertical heat flux
end

subsection Particles
set Data output format = ascii, vtu
set Data output format = vtu
set List of particle properties = position, velocity, initial composition, pT path, initial position
set Number of particles = 100000
set Number of grouped files = 1
set Update ghost particles = true
set Particle generator name = uniform radial
set Particle generator name = reference cell
set Time between data output = 1e6
set Allow cells without particles = true
set Interpolation scheme = cell average
set Write in background thread = true

set Load balancing strategy = remove and add particles
set Minimum particles per cell = 5
set Minimum particles per cell = 80

subsection Generator
subsection Uniform radial
set Minimum radius = 3480000
set Maximum radius = 6370000
set Radial layers = 50
subsection Reference cell
set Number of particles per cell per direction = 3
end
end

end
end

0 comments on commit 9814c9f

Please sign in to comment.