Subduction 3D miniapp

Project.toml

@@ -31,7 +31,7 @@ AMDGPU = "0.6, 0.7, 0.8"
 Adapt = "3"
 CUDA = "4.4.1, 5"
 CellArrays = "0.2"
-GeoParams = "0.5.8"
+GeoParams = "0.6"
 HDF5 = "0.17.1"
 ImplicitGlobalGrid = "0.15.0"
 MPI = "0.20"

Subduction_GMG_2D.jl

@@ -0,0 +1,58 @@
+#=
+# Creating 2D numerical model setups
+
+### Aim
+The aim of this tutorial is to show you how to create 2D numerical model setups that can be used as initial setups for other codes.
+
+=#
+
+
+#=
+### 2D Subduction setup
+
+Lets start with creating a 2D model setup in cartesian coordinates, which uses the `CartData` data structure
+=#
+using GeophysicalModelGenerator
+
+function GMG_subduction_2D(nx, ny)
+    # Our starting basis is the example above with ridge and overriding slab
+    nx, nz = nx, ny
+    x      = range(-1000,1000, nx);
+    z      = range(-660,0,    nz);
+    Grid2D = CartData(xyz_grid(x,0,z))
+    Phases = zeros(Int64, nx, 1, nz);
+    Temp   = fill(1350.0, nx, 1, nz);
+    lith   = LithosphericPhases(Layers=[15 20 55], Phases=[3 4 5], Tlab=1250)
+    # mantle = LithosphericPhases(Phases=[1])
+
+    # add_box!(Phases, Temp, Grid2D; xlim=(-1000, 1000), zlim=(-600.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
+    # Phases .= 0
+
+    # Lets start with defining the horizontal part of the overriding plate. Note that we define this twice with different thickness to deal with the bending subduction area:
+    add_box!(Phases, Temp, Grid2D; xlim=(200,1000), zlim=(-150.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
+    add_box!(Phases, Temp, Grid2D; xlim=(0,200), zlim=(-50.0, 0.0), phase = lith, T=HalfspaceCoolingTemp(Age=80));
+
+    # The horizontal part of the oceanic plate is as before:
+    v_spread_cm_yr = 3      #spreading velocity
+    lith = LithosphericPhases(Layers=[15 55], Phases=[1 2], Tlab=1250)
+    add_box!(Phases, Temp, Grid2D; xlim=(-1000,0.0), zlim=(-150.0, 0.0), phase = lith, T=SpreadingRateTemp(SpreadingVel=v_spread_cm_yr));
+
+    # Yet, now we add a trench as well. The starting thermal age at the trench is that of the horizontal part of the oceanic plate:
+    AgeTrench_Myrs = 1000e3/(v_spread_cm_yr/1e2)/1e6    #plate age @ trench
+
+    # We want to add a smooth transition from a halfspace cooling 1D thermal profile to a slab that is heated by the surrounding mantle below a decoupling depth `d_decoupling`.
+    T_slab = LinearWeightedTemperature( F1=HalfspaceCoolingTemp(Age=AgeTrench_Myrs), F2=McKenzie_subducting_slab(Tsurface=0,v_cm_yr=v_spread_cm_yr, Adiabat = 0.0))
+
+    # in this case, we have a more reasonable slab thickness: 
+    trench = Trench(Start=(0.0,-100.0), End=(0.0,100.0), Thickness=100.0, θ_max=30.0, Length=600, Lb=200, 
+                    WeakzoneThickness=15, WeakzonePhase=6, d_decoupling=125);
+    add_slab!(Phases, Temp, Grid2D, trench, phase = lith, T=T_slab);
+
+    # Lithosphere-asthenosphere boundary:
+    ind = findall(Temp .> 1250 .&& (Phases.==2 .|| Phases.==5));
+    Phases[ind] .= 0;
+
+    Grid2D = addfield(Grid2D,(;Phases, Temp))
+
+end
+

_typos.toml

@@ -2,3 +2,4 @@
 Ths = "Ths"
 oly = "oly"
 iy = "iy"
+nd = "nd"

docs/src/man/ShearBands.md

@@ -10,7 +10,6 @@ using JustRelax, JustRelax.JustRelax2D, JustRelax.DataIO
 const backend_JR = CPUBackend
 ```
 
-
 We will also use `ParallelStencil.jl` to write some device-agnostic helper functions:
 ```julia
 using ParallelStencil
@@ -54,7 +53,7 @@ dt      = η0/G0/4.0     # assumes Maxwell time of 4
 el_bg   = ConstantElasticity(; G=G0, Kb=4)
 el_inc  = ConstantElasticity(; G=Gi, Kb=4)
 visc    = LinearViscous(; η=η0)
-pl      = DruckerPrager_regularised(;  # non-regularized plasticity
+pl      = DruckerPrager_regularised(;  # regularized plasticity
     C    = C,
     ϕ    = ϕ,
     η_vp = η_reg,

docs/src/man/equations.md

@@ -1,7 +1,55 @@
+# Stokes equations
+
+Stokes equations for compressible flow are described by
+
+$$\nabla\cdot\boldsymbol{\tau} + \nabla p = \boldsymbol{f} $$
+
+$$\nabla\cdot\boldsymbol{v} + \beta \frac{\partial p}{\partial t} + \alpha \frac{\partial T}{\partial t} = 0$$
+
+where $\nabla = \left(\frac{\partial }{\partial x_i}, ..., \frac{\partial }{\partial x_n} \right)$ is the nabla operator, $\boldsymbol{\tau}$ is the deviatoric stress tensor, $p$ is the pressure, $\boldsymbol{f}$ is the body forces vector, $\boldsymbol{v}$ is the velocity field, and $\beta$ is in the inverse of the bulk modulus. In the simple case of a linear rheology, the constitutive equation for an isotropic flow is $\boldsymbol{\tau} = 2\eta\dot{\boldsymbol{\varepsilon}}$, where $\dot{\boldsymbol{\varepsilon}}$ is the deviatoric strain rate tensor. Heat diffusion
+
+$$\rho C_p \frac{\partial T}{\partial t} = - \nabla q + Q$$
+
+$$q = - K \nabla T$$
+
+where $\rho$ is the density, $C_p$ is the heat diffusion, $K$ is the heat conductivity, $Q$ is the sum of any amount of source terms, and $T$ is the temperature.
+
 # Pseudo-transient iterative method
 
+The pseudo-transient method consists in augmenting the right-hand-side of the target PDE with a pseudo-time derivative (where $\psi$ is the pseudo-time) of the primary variables. 
+
+## Stokes
+
+The pseudo-transient formulation of the Stokes equations yields:
+
+$$\widetilde{\rho}\frac{\partial \boldsymbol{u}}{\partial \psi} + \nabla\cdot\boldsymbol{\tau} - \nabla p = \boldsymbol{f}$$
+
+$$\frac{1}{\widetilde{K}}\frac{\partial p}{\partial \psi} + \nabla\cdot\boldsymbol{v} = \beta \frac{\partial p}{\partial t} + \alpha \frac{\partial T}{\partial t}$$
+
+We also do a continuation on the constitutive law:
+
+$$\frac{1}{2\widetilde{G}} \frac{\partial\boldsymbol{\tau}}{\partial\psi}+ \frac{1}{2G}\frac{D\boldsymbol{\tau}}{Dt} + \frac{\boldsymbol{\tau}}{2\eta} = \dot{\boldsymbol{\varepsilon}}$$
+
+where the wide tile denotes the effective damping coefficients and $\psi$ is the pseudo-time step. These are defined as in [Raess et al, 2022](https://doi.org/10.5194/gmd-2021-411):
+
+$$\widetilde{\rho} = Re\frac{\eta}{\widetilde{V}L}, \qquad
+  \widetilde{G} = \frac{\widetilde{\rho} \widetilde{V}^2}{r+2}, \qquad
+  \widetilde{K} = r \widetilde{G}$$
+
+and
+
+$$\widetilde{V} = \sqrt{ \frac{\widetilde{K} +2\widetilde{G}}{\widetilde{\rho}}}, \qquad
+    r = \frac{\widetilde{K}}{\widetilde{G}}, \qquad
+    Re = \frac{\widetilde{\rho}\widetilde{V}L}{\eta}$$
+
+where the P-wave $\widetilde{V}=V_p$ is the characteristic velocity scale for Stokes, and $Re$ is the Reynolds number. 
+
 ## Heat diffusion
 
-## Stokes equations
+The pseudo-transient heat-diffusion equation is:
+
+$$\widetilde{\rho}\frac{\partial T}{\partial \psi} + \rho C_p \frac{\partial T}{\partial t} = \nabla(K\nabla T) = -\nabla q$$
+
+We use a second order pseudo-transient scheme were continuation is also done on the flux, so that:
 
-## Constitutive equations
+$$\widetilde{\theta}\frac{\partial q}{\partial \psi} + q  = -K\nabla T$$

docs/src/man/subduction2d/rheology.md

@@ -0,0 +1,76 @@
+# Using GeoParams.jl to define the rheology of the material phases
+
+We will use the same physical parameters as the ones defined in [Hummel, et al 2024](https://doi.org/10.5194/se-15-567-2024).
+
+The thermal expansion coefficient $\alpha$ and heat capacity $C_p$ are the same for all the material phases
+
+```julia
+α  = 2.4e-5 # 1 / K
+Cp = 750    # J / kg K
+```
+
+The density of all the phases is constant, except for the oceanic lithosphere, which uses the pressure and temperature dependent equation of state $\rho = \rho_0 \left(1 - \alpha (T-T_0) - \beta (P-P_0) \right)$, where $\rho_0 = \rho (T=1475 \text{C}^{\circ})=3200 \text{kg/m}^3$.which corresponds to the `PT_Density` object from `GeoParams.jl`:
+
+```julia
+density_lithosphere = PT_Density(; ρ0=3.2e3, α = α, β = 0e0, T0 = 273+1474)
+```
+
+We will run the case where the rheology is given by a combination of dislocation and diffusion creep for wet olivine, 
+
+```julia
+disl_wet_olivine  = SetDislocationCreep(Dislocation.wet_olivine1_Hirth_2003)
+diff_wet_olivine  = SetDiffusionCreep(Diffusion.wet_olivine_Hirth_2003)
+```
+
+and where plastic failure is given by the formulation from [Duretz et al, 2021](https://doi.org/10.1029/2021GC009675)
+```julia
+# non-regularized plasticity
+ϕ               = asind(0.1)
+C               = 1e6        # Pa
+plastic_model   = DruckerPrager_regularised(; C = C, ϕ = ϕ, η_vp=η_reg, Ψ=0.0) 
+```    
+
+Finally we define the rheology objects of `GeoParams.jl`
+```julia
+rheology = (
+    SetMaterialParams(;
+        Name              = "Asthenoshpere",
+        Phase             = 1,
+        Density           = ConstantDensity(; ρ=3.2e3),
+        HeatCapacity      = ConstantHeatCapacity(; Cp = Cp),
+        Conductivity      = ConstantConductivity(; k  = 2.5),
+        CompositeRheology = CompositeRheology( (LinearViscous(; η=1e20),)),
+        Gravity           = ConstantGravity(; g=9.81),
+    ),
+    SetMaterialParams(;
+        Name              = "Oceanic lithosphere",
+        Phase             = 2,
+        Density           = density_lithosphere,
+        HeatCapacity      = ConstantHeatCapacity(; Cp = Cp),
+        Conductivity      = ConstantConductivity(; k  = 2.5),
+        CompositeRheology = CompositeRheology( 
+            (
+                disl_wet_olivine,
+                diff_wet_olivine,
+                plastic_model,
+            ) 
+        ),
+    ),
+    SetMaterialParams(;
+        Name              = "oceanic crust",
+        Phase             = 3,
+        Density           = ConstantDensity(; ρ=3.2e3),
+        HeatCapacity      = ConstantHeatCapacity(; Cp = Cp),
+        Conductivity      = ConstantConductivity(; k  = 2.5),
+        CompositeRheology = CompositeRheology( (LinearViscous(; η=1e20),)),
+    ),
+    SetMaterialParams(;
+        Name              = "StickyAir",
+        Phase             = 4,
+        Density           = ConstantDensity(; ρ=1), # water density
+        HeatCapacity      = ConstantHeatCapacity(; Cp = 1e34),
+        Conductivity      = ConstantConductivity(; k  = 3),
+        CompositeRheology = CompositeRheology((LinearViscous(; η=1e19),)),
+    ),
+)
+```

docs/src/man/subduction2d/setup.md

@@ -0,0 +1,101 @@
+# Model setup
+As described in the original [paper](https://doi.org/10.5194/se-15-567-2024), the domain consists of a Cartesian box of $\Omega \in [0, 3000] \times [0, -660]$ km, with two 80km thick oceanic plates over the asthenospheric mantle. 
+
+We will use GeophysicalModelGenerator.jl to generate the initial geometry, material phases, and thermal field of our models. We will start by defining the dimensions and resolution of our model, as well as initializing the `Grid2D` object and two arrays `Phases` and `Temp` that host the material phase (given by an integer) and the thermal field, respectively.
+
+```julia
+nx, nz        = 512, 218 # number of cells per dimension
+Tbot          = 1474.0   # [Celsius]
+model_depth   = 660      # [km]
+air_thickness = 10       # [km]
+Lx            = 3000     # model length [km]
+x             = range(0, Lx, nx);
+z             = range(-model_depth, air_thickness, nz);
+Grid2D        = CartData(xyz_grid(x,0,z))
+Phases        = zeros(Int64, nx, 1, nz);
+Temp          = fill(Tbot, nx, 1, nz);
+```
+
+In this model we have four material phases with their respective phase numbers: 
+
+| Material            | Phase number |
+| :----------------   | :----------: |
+| asthenosphere       |       0      |
+| oceanic lithosphere |       1      |
+| oceanic crust       |       3      |
+| sticky air          |       4      |
+
+We will start by initializing the model as asthenospheric mantle, with a thermal profile given by the half-space cooling model with an age of 80 Myrs.
+
+```julia
+add_box!(
+    Phases, 
+    Temp, 
+    Grid2D; 
+    xlim    = (0, Lx),
+    zlim    = (-model_depth, 0.0), 
+    phase   = LithosphericPhases(Layers=[], Phases=[0]), 
+    T       = HalfspaceCoolingTemp(Tsurface=20, Tmantle=Tbot, Age=80,Adiabat=0.4)
+)
+```
+![](setup_1.png)
+
+Next we add a horizontal 80km thick oceanic lithosphere. Note that we leave a 100km buffer zone next to the vertical boundaries of the domain, to facilitate the sliding of the oceanic plates.
+```julia
+add_box!(
+    Phases, 
+    Temp, 
+    Grid2D; 
+    xlim    = (100, Lx-100), # 100 km buffer zones on both sides
+    zlim    = (-model_depth, 0.0),
+    phase   = LithosphericPhases(Layers=[80], Phases=[1 0]), 
+    T       = HalfspaceCoolingTemp(Tsurface=20, Tmantle=Tbot, Age=80, Adiabat=0.4)
+)
+```
+![](setup_2.png)
+
+As in the original paper, we add a 8km thick crust on top of the subducting oceanic plate.
+```julia
+# Add right oceanic plate crust
+add_box!(
+    Phases, 
+    Temp, 
+    Grid2D; 
+    xlim    = (Lx-1430, Lx-200), 
+    zlim    = (-model_depth, 0.0), 
+    Origin  = nothing, StrikeAngle=0, DipAngle=0,
+    phase   = LithosphericPhases(Layers=[8 72], Phases=[2 1 0]), 
+    T       = HalfspaceCoolingTemp(Tsurface=20, Tmantle=Tbot, Age=80, Adiabat=0.4)
+)
+```
+![](setup_3.png)
+
+And finally we add the subducting slab, with the trench located at 1430km from the right-hand-side boundary.
+
+```julia
+add_box!(
+    Phases, 
+    Temp, 
+    Grid2D; 
+    xlim    = (Lx-1430, Lx-1430-250), 
+    zlim    = (-80, 0.0), 
+    Origin  = (nothing, StrikeAngle=0, DipAngle=-30),
+    phase   = LithosphericPhases(Layers=[8 72], Phases=[2 1 0]), 
+    T       = HalfspaceCoolingTemp(Tsurface=20, Tmantle=Tbot, Age=80, Adiabat=0.4)
+)
+```
+![](setup_4.png)
+
+```julia
+surf = Grid2D.z.val .> 0.0 
+@views Temp[surf] .= 20.0
+@views Phases[surf] .= 3
+```
+![](setup_5.png)
+
+```julia
+li     = (abs(last(x)-first(x)), abs(last(z)-first(z))) .* 1e3 # in meters
+origin = (x[1], z[1]) .* 1e3 # lower-left corner of the domain
+Phases = Phases[:,1,:] .+ 1  # +1 because Julia is 1-indexed
+Temp   = Temp[:,1,:].+273    # in Kelvin
+```