contribute cookbook: simple subduction initiation

cookbooks/subduction_initiation/doc/subduction_initiation.md

@@ -0,0 +1,91 @@
+(sec:cookbooks:subd-init)=
+# Subduction initiation from Matsumoto and Tomoda (1983)
+
+*This section was contributed by Cedric Thieulot.*
+
+The setup for this experiment originates in {cite:t}`matsumoto:tomoda:1983`.
+
+In this very early computational geodynamics paper the authors
+are interested in predicting the future process of crustal and lithospheric movement
+at the gigantic fracture zone in the Northeastern
+Pacific by means of numerical simulation
+of contact between two kinds of viscous fluid of different density.
+They then solve the incompressible isothermal linear Stokes equations
+with a stream function approach and the resulting equation is solved by means of the
+Finite Difference method and the Marker and Cell (MAC) method.
+
+The domain is a 2D Cartesian box of size $L_x \times L_y=(400~\text{ km},180~\text{ km})$.
+There are three fluids in the domain: water ($\rho_w=1030~\text{ kg m}^{-3}$, $\eta_w=10^{-3}~\text{ Pa s}$),
+lithosphere ($\rho_l=3300~\text{ kg m}^{-3}$, $\eta_l$) and
+asthenosphere ($\rho_a=3200~\text{ kg m}^{-3}$, $\eta_a$).
+Note that although the water viscosity is correct, this is probably not
+the value that was used in the authors since it would most likely lead to
+numerical errors. We then set $\eta_w=10^{19}~\text{ Pa s}$ (we can then speak of 'sticky water').
+Also, all viscosities in the paper
+are expressed in Poise with $1~\text{ Poise}=0.1~\text{ Pa s}$.
+
+The setup is shown in {numref}`fig:subduction-initiation-setup`
+and the list of simulations run by the authors is in
+{numref}`tab:quickref`.
+
+```{figure-md} fig:subduction-initiation-setup
+<img src="setup.*" width="70%" />
+
+$D_w= 10~\text{ km}$,
+$D_l= 50~\text{ km}$,
+$d_w=8~\text{ km}$,
+$d_l=10~\text{ km}$. These values
+represent the easternmost part of the Mendocino
+Fracture Zone.
+Taken from {cite:t}`matsumoto:tomoda:1983`.
+```
+
+```{table} List of all cases
+:name: tab:quickref
+| Case                 | $\eta_l~(\text{ Pa s})$      | $\eta_a~(\text{ Pa s})$       | domain size  (km)|
+| :------------------- | :-----------: | :------------: | :---------: |
+1 | $10^{22}$ | $10^{21}$ | $400\times 180$ |
+2 | $10^{22}$ | $10^{20}$ | $400\times 180$ |
+3 | $10^{22}$ | $10^{19}$ | $400\times 180$ |
+4 | $10^{23}$ | $10^{21}$ | $400\times 180$ |
+5 | $10^{22}$ | $10^{20}$ | $800\times 140$ |
+6 | $10^{22}$ | $10^{19}$ | $800\times 140$ |
+```
+
+
+Rather interestingly the model is built in such a way that the
+lithostatic pressure is uniform at the bottom of the domain.
+Models are run for 50 Myr.
+Boundary conditions are free-slip on all sides of the domain.
+Unfortunately the authors do not specify the mesh resolution that was used
+but {numref}`fig:subduction-initiation-results1` gives us an idea.
+
+```{figure-md} fig:subduction-initiation-results1
+<img src="results1.*" width="90%" />
+
+Computer output of the result of Case 6(c). \# indicates a particle
+representing the material of older lithosphere, * is that of
+younger lithosphere, : is older asthenosphere and $\cdot$ is younger asthenosphere.
+Taken from {cite:t}`matsumoto:tomoda:1983`.
+```
+
+Two parameter files are provided for Case 1: one using
+compositional fields, and one using the particle-in-cell approach.
+Other cases can be run by changing the viscosities and/or the domain size accordingly in the
+parameter files.
+Results for Case 1 in the original publication are shown in {numref}`fig:subduction-initiation-results2`
+and results obtained in ASPECT are shown in {numref}`fig:subduction-initiation-results3`.
+
+```{figure-md} fig:subduction-initiation-results2
+<img src="results2.*" width="70%" />
+
+Results of calculation in Case 1. (a) 0 Ma; (b) 12 Ma; (c) 21.4 Ma; (d) 30.7 Ma; (e) 46.7 Ma.
+Taken from {cite:t}`matsumoto:tomoda:1983`.
+```
+
+```{figure-md} fig:subduction-initiation-results3
+<img src="results3.*" width="95%" />
+
+Case 1 (obtained with the compositional fields approach) at
+(a) 0 Ma; (b) 12 Ma; (c) 21 Ma; (d) 31 Ma; (e) 41 Ma.
+```

cookbooks/subduction_initiation/subduction_initiation_compositional_fields.prm

@@ -0,0 +1,98 @@
+#parameter file for replicating Matsumoto & Tomoda 1983
+
+set Dimension                              = 2
+set Start time                             = 0
+set End time                               = 50e6
+set Use years in output instead of seconds = true
+set CFL number                             = 0.25
+set Output directory                       = output-subduction-initiation-comp
+set Pressure normalization                 = surface
+
+subsection Solver parameters
+  subsection Stokes solver parameters
+    set Number of cheap Stokes solver steps = 0
+  end
+end
+
+subsection Geometry model
+  set Model name = box
+  subsection Box
+    set X extent  = 400e3
+    set Y extent  = 180e3
+    set X repetitions = 2
+  end
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = left, right, bottom, top
+end
+
+
+#materials are:
+# asthenosphere left,
+# asthenosphere right,
+# lithosphere left,
+# lithosphere right,
+# water.
+
+subsection Material model
+  set Model name = multicomponent
+  subsection Multicomponent
+    set Densities = 3200, 3200, 3300, 3300, 1030
+    set Viscosities = 1e21, 1e21, 1e22, 1e22, 1e19
+    set Viscosity averaging scheme = harmonic
+    set Thermal expansivities = 0
+  end
+end
+
+subsection Gravity model
+  set Model name = vertical
+  subsection Vertical
+    set Magnitude = 9.81
+  end
+end
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators = bottom, top
+  set List of model names = box
+end
+
+subsection Initial temperature model
+  set Model name = function
+  subsection Function
+    set Function expression = 0
+  end
+end
+
+subsection Compositional fields
+  set Number of fields = 4
+end
+
+subsection Initial composition model
+  set Model name = function
+  subsection Function
+    set Variable names      = x,y
+    set Function constants  = L0=300e3
+    set Function expression = if((x>=L0 && y<162e3),1,0) ;\
+                              if((x<=L0 && y>120e3 && y<=170e3),1,0) ;\
+                              if((x>=L0 && y>=162e3 && y<=172e3),1,0) ;\
+                              if((x<=L0 && y>170e3) || (x>=L0 && y>172e3),1,0)
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement   = 2
+  set Initial global refinement     = 5
+  set Refinement fraction           = 0.9
+  set Strategy                      = composition
+  set Coarsening fraction           = 0
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization, velocity statistics, composition statistics, pressure statistics, material statistics, global statistics
+  subsection Visualization
+    set List of output variables = density, viscosity, strain rate
+    set Time between graphical output = 0
+    set Interpolate output = false
+  end
+end

cookbooks/subduction_initiation/subduction_initiation_particle_in_cell.prm

@@ -0,0 +1,104 @@
+#parameter file for replicating Matsumoto & Tomoda 1983
+
+set Dimension                              = 2
+set Start time                             = 0
+set End time                               = 50e6
+set Use years in output instead of seconds = true
+set CFL number                             = 0.25
+set Output directory                       = output-subduction-initiation-pic
+set Pressure normalization                 = surface
+
+#subsection Solver parameters
+#  subsection Stokes solver parameters
+#    set Number of cheap Stokes solver steps = 0
+#  end
+#end
+
+subsection Geometry model
+  set Model name = box
+  subsection Box
+    set X extent  = 400e3
+    set Y extent  = 180e3
+    set X repetitions = 2
+  end
+end
+
+subsection Boundary velocity model
+  set Tangential velocity boundary indicators = left, right, bottom, top
+end
+
+subsection Material model
+  set Model name = multicomponent
+  subsection Multicomponent
+    set Densities = 0, 3200, 3200, 3300, 3300, 1030
+    set Viscosities = 1e30, 1e21, 1e21, 1e22, 1e22, 1e19
+    set Viscosity averaging scheme = harmonic
+    set Thermal expansivities = 0
+  end
+  set Material averaging = harmonic average
+end
+
+subsection Gravity model
+  set Model name = vertical
+  subsection Vertical
+    set Magnitude = 9.81
+  end
+end
+
+subsection Boundary temperature model
+  set Fixed temperature boundary indicators = bottom, top
+  set List of model names = box
+end
+
+subsection Initial temperature model
+  set Model name = function
+  subsection Function
+    set Function expression = 0
+  end
+end
+
+subsection Compositional fields
+  set Number of fields = 5
+  set Names of fields = asth_left, asth_right, left_lith, right_lith, water
+  set Compositional field methods = particles, particles, particles, particles, particles
+  set Mapped particle properties = asth_left:initial asth_left, asth_right:initial asth_right, left_lith:initial left_lith, right_lith:initial right_lith, water: initial water
+end
+
+subsection Initial composition model
+  set Model name = function
+  subsection Function
+    set Variable names      = x,y
+    set Function constants  = L0=300e3
+    set Function expression = if((x<=L0 && y<=120e3),1,0) ;\
+                              if((x>=L0 && y<162e3),1,0) ;\
+                              if((x<=L0 && y>120e3 && y<=170e3),1,0) ;\
+                              if((x>=L0 && y>=162e3 && y<=172e3),1,0) ;\
+                              if((x<=L0 && y>170e3) || (x>=L0 && y>172e3),1,0)
+  end
+end
+
+subsection Mesh refinement
+  set Initial adaptive refinement   = 0
+  set Initial global refinement     = 6
+  set Refinement fraction           = 0.9
+  set Strategy                      = composition
+  set Coarsening fraction           = 0
+end
+
+subsection Postprocess
+  set List of postprocessors = visualization, velocity statistics, composition statistics, pressure statistics, material statistics, global statistics, particles
+  subsection Particles
+    set Number of particles      = 350000
+    set Time between data output = 0
+    set Data output format       = vtu
+    set List of particle properties = initial composition, velocity
+    set Particle generator name = random uniform
+    set Interpolation scheme = cell average #default
+    set Update ghost particles = true
+  end
+  subsection Visualization
+    set List of output variables = density, viscosity, strain rate
+    set Time between graphical output = 0
+    set Interpolate output = false
+  end
+end

doc/modules/changes/20230712_thieulot

@@ -0,0 +1,5 @@
+New: There is now a new cookbook which is a very 
+simple subduction initiation model as published 
+by Matsumoto and Tomoda (1983). 
+<br>
+(Cedric Thieulot, 2023/07/12)

doc/sphinx/references.bib

@@ -22,6 +22,16 @@ @article{halter:macherel:schmalholz:2022
   year={2022}
 }
 
+@article{matsumoto:tomoda:1983,
+  title={Numerical simulation of the initiation of subduction at the fracture zone},
+  author={Matsumoto, Takeshi and Tomoda, Yoshibumi},
+  journal={Journal of Physics of the Earth},
+  volume={31},
+  number={3},
+  pages={183--194},
+  doi={10.4294/jpe1952.31.183},
+  year={1983}
+}
 
 @article{kiefer:hager:1992,
   title={{Geoid anomalies and dynamic topography from convection in cylindrical geometry: applications to mantle plumes on Earth and Venus}},

doc/sphinx/user/cookbooks/geophysical-setups.md

@@ -106,5 +106,6 @@ cookbooks/global_melt/doc/global_melt.md
 cookbooks/mid_ocean_ridge/doc/mid_ocean_ridge.md
 cookbooks/kinematically_driven_subduction_2d/doc/kinematically_driven_subduction_2d.md
 cookbooks/inclusions/doc/inclusions.md
+cookbooks/subduction_initiation/doc/subduction_initiation.md
 cookbooks/future/README.md
 :::