Tags: category:cookbook feature:3d feature:cartesian feature:particles feature:compositional-fields

Olivine Fabric Developments Under Simple Shear

This section was contributed by Xiaochuan Tian, with help from Yijun Wang and Menno Fraters. It is based on a section in Fraters and Billen [2021] by Menno Fraters and Magali Billen published in 2021.

This cookbook explains how to set up a numerical experiment for fabric developments of a single olivine particle under simple shear macroscopic strain. It uses ASPECT crystal preferred orientation (CPO) implementation, which is described in detail in Fraters and Billen [2021]. The fabric calculation is based on DRex Kaminski et al. [2004]. This cookbook describes how to use ASPECT’s CPO implementation to calculate the major olivine fabrics observed in the lab experiments. It focuses on reproducing Fig. 3 of [Fraters and Billen, 2021].

Motivation

The Earth’s plates can move relative to each other because the underlying mantle can be deformed to accommodate such motions. This motion reflects dynamics inside the planet, which cannot be observed directly, as the depth of interest is beyond human drilling capability. However, indirect observations from seismic anisotropy provide information about the directionality of seismic wave speeds propagating inside the Earth. If such directionality is controlled by the mantle’s motion or deformation, we can use seismic anisotropy observations to infer mantle deformation.

It is important to note that this cookbook does not directly compare model outputs with seismic anisotropy observations. Instead, it provides the full elastic tensor, including its symmetry axes and fast directions, at the locations of particles in the model. To make direct comparisons with seismic observations, you would still need to process these results through a separate code that simulates seismic wave propagation through the modeled medium.

Olivine is a major mineral of the Earth’s upper mantle where continuous deformation takes place. People care about olivine fabric developments under simple shear as it may provide a bridge that links seismic anisotropy observations to flow in the upper mantle that accommodates plate motions. Indeed, high temperature and pressure lab experiments using Griggs apparatus to investigate olivine fabrics developments under simple shear do find systematic types of fabrics based on the stress and water content conditions (see Karato et al. [2008] for details). Note that the effects of secondary phases, such as enstatite Kaminski et al. [2004] can also be included in the method, we here only focus on olivine for simplicity. Meanwhile, following the example input file from Fraters and Billen [2021],the input file provided here also includes an example of enstatite calculation with a volume fraction of 0.3.

Model setup

Following Fraters and Billen [2021], we prescribe simple shear in a 3d Cartesian box/cube with dimensions of \(1 \times 1 \times 1 \) \([m^3]\). The shear strain rate is set to \(\dot{\epsilon}_{xz} = -5\times 10^{-6} [s^{-1}]\). The particle with olivine/enstatite volume fraction of 0.7/0.3 is placed right at the center of the cubic box so it stays stationary. The DRex implementation keep tracks of rotations of crystal grains within the particle under macroscopic deformation.

Fig. 178 3D shear box with velocity vectors. The grey ball is the central olivine particle.

The model computes how crystal grains rotate and align under simple shear and the pole figures visualize this rotation.

The input file

Rather than solving the Stokes equation with boundary conditions, we instead prescribe a constant simple shear strain rate field by setting “Nonlinear solver scheme” to “single Advection, no Stokes” and prescribing the Stokes solution with a function. In this case, the only model property that we care about is the prescribed strain rate that is exerted onto the olivine particle. Parameters useful for the Stokes equation is not important in this case.

# we prescribe the simple shear strain rate by setting solver scheme
# as "single Advection, no Stokes"
# see below "subsection Prescribed Stokes solution "
set Nonlinear solver scheme  = single Advection, no Stokes
# The model run to 3e5 seconds, at which the total strain will
# be 1.5 with strain rate of 5e-6 1/s.
set End time                 = 3e5
set Use years instead of seconds           = false

# Set zero gravity
subsection Gravity model
  set Model name = vertical

  subsection Vertical
    set Magnitude = 0
  end
end

# 1m by 1m by 1m cube
# center is at (0,0,0)
subsection Geometry model
  set Model name = box

  subsection Box
    set X extent = 1
    set Y extent = 1
    set Z extent = 1
    set Box origin X coordinate = -0.5
    set Box origin Y coordinate = -0.5
    set Box origin Z coordinate = -0.5
  end
end

# Here prescribing velocity field so that it results in
# simple shear strain rate epsilon_xz = -5e6 1/s
# as z ranges from -0.5 to 0.5 meter.
subsection Prescribed Stokes solution
  set Model name = function

  subsection Velocity function
    set Variable names = x,y,z,t
    set Function expression = z/1e5;0;0 ## Annotation set velocity condition
  end
end

# The default X,Y,Z repetitions is 1, so will be 2 by 2 by 2 grids without refinement
# Here with setting "Initial global refinement" to be 1, it becomes
# 4 by 4 by 4 grids.
subsection Mesh refinement
  set Initial global refinement                = 1
  set Initial adaptive refinement              = 0
  set Time steps between mesh refinement       = 10000
end

The second part of the input file specifies what to output and for what minerals to keep track of CPO development. It also specifies the parameters for the DRex algorithm.

# Specify postprocessor objects that will be run at the end of each time step
# The particles and crystal preferred orientation are required for CPO calculations
subsection Postprocess
  set List of postprocessors = velocity statistics, visualization, particles, crystal preferred orientation

  subsection Visualization
    set Time between graphical output = 1e5
    set List of output variables = stress, material properties, strain rate
  end

  # Specify CPO data output
  subsection Crystal Preferred Orientation
    set Time between data output = 1e5
    set Write in background thread = true
    set Compress cpo data files = false # so results are human readable
    set Write out raw cpo data = mineral 0: volume fraction, mineral 0: Euler angles #, mineral 1: volume fraction, mineral 1: Euler angles #Specify what to output, here will only output volume fraction and Euler angles for mineral 0 (olivine)
     # Specify output parameters for draw volume weighted cpo, the pole figures are plotted with these data
    set Write out draw volume weighted cpo data = mineral 0: volume fraction, mineral 0: Euler angles, mineral 1: volume fraction, mineral 1: Euler angles
  end

  subsection Particles
    set Time between data output = 1e4
    set Data output format       = vtu
    # Exclude output properties so as to allow fast reading by paraview
    # because if number of particles is large, paraview will crash
    set Exclude output properties = volume fraction, rotation_matrix
  end
end

subsection Particles
  set List of particle properties = integrated strain, crystal preferred orientation, cpo elastic tensor, cpo bingham average, integrated strain invariant, velocity, pT path

  # use an ASCII file particle_one.dat to specify initial locations of particles
  set Particle generator name = ascii file

  subsection Generator
    subsection Ascii file
      # you may need to change the below path to the absolute path in your local computer
      # the particle_one.dat is located in the cookbook folder
      set Data directory = ./
      set Data file name = particle_one.dat
    end
  end

  # parameters for DRex
  subsection Crystal Preferred Orientation
    # a random seed to 'randomize' initial Euler angles for the grains
    set Random number seed = 301
    set Number of grains per particle = 500
    set Property advection method = Backward Euler
    set Property advection tolerance = 1e-15
    set CPO derivatives algorithm = D-Rex 2004

    subsection Initial grains
      # prescribe types of mineral
      # The A,B,C,D,E types of olivine fabric can be set here by changing
      # A-fabric to B-fabric for example
      set Minerals = Olivine: A-fabric, Enstatite

      # The two minerals are not interacting with each other.
      set Volume fractions minerals = 0.7,0.3
    end

    # DRex parameters following Kaminski et al., 2004
    subsection D-Rex 2004
      set Mobility = 125
      set Stress exponents = 3.5
      set Exponents p = 1.5
      set Nucleation efficiency = 5
      set Threshold GBS = 0.3
    end
  end

  subsection CPO Bingham Average
    set Random number seed = 200
  end
end

The complete input file is located at /cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/olivineA.prm.

Plotting pole figures

To plot pole figures from the output in out_drex_olivineA/particles_cpo, you could use MTEX or cpo_analyzer. For example, inside the cpo_analyzer folder, run:

$ cargo run --release path/to/cpo_pole_figure_config.toml

A folder named ‘CPO_figures’ will be generated in the folder ‘out_drex_olivineA’ that contains the pole figure results. An example config file can be found here: /cookbooks/crystal_preferred_orientation_olivine_fraters_billen_2021/cpo_pole_figure_config.toml.

Model results

Fig. 179 Initial randomized CPO at time zero.

By changing “Olivine: A-fabric” to “Olivine: B-fabric” you can generate olivine fabric results for B-type and other types as well shown below.

Fig. 180 Olivine A-E type fabrics under simple shear with a shear strain of 1.5 at 3e5 seconds of model time under constant shear strain rate: \(\dot{\epsilon}_{xz} = -5\times 10^{-6} [s^{-1}]\). Note that the color scales range from 0 to 5 but actual results may be larger.

When we look at the visualization output of this model (see also Fig. 180), we can see the results reproduce the pole figures in Figure 3 of Fraters and Billen [2021]. There is a minor difference in the type D fabric (the girdle is less pronounced) when compared to Fig. 3 of (Karato et al., 2008), which could be investigated in more details in the future.

Extending the model

This model of a unit cube under simple shear can be extended in various ways. Some ideas for adapting the model setup are:

• Setup high pressure and temperature lab experiments for olivine CPO. After benchmarking against the lab experiment results, it could then be extended beyond the parameter space limited by the lab experiment.

• Add CPO calculations in large scale geodynamic models like mantle plumes, lower crust delamination and mid-ocean ridges.

• Include anisotropic viscosity to feed this CPO back to influence the stokes solution.

References:

• Fraters, M. R. T., & Billen, M. I. (2021). On the implementation and usability of crystal preferred orientation evolution in geodynamic modeling. Geochemistry, Geophysics, Geosystems, 22(10), e2021GC009846.

• Kaminski, E., Ribe, N. M., & Browaeys, J. T. (2004). D-Rex, a program for calculation of seismic anisotropy due to crystal lattice preferred orientation in the convective upper mantle. Geophysical Journal International, 158(2), 744–752.

• Karato, S., Jung, H., Katayama, I., & Skemer, P. (2008). Geodynamic significance of seismic anisotropy of the upper mantle: New insights from laboratory studies. Annu. Rev. Earth Planet. Sci., 36, 59–95.