Unphysical viscosity structure caused by initial compositional distribution

[XHL-CUG]

Hello everyone,

I have encountered an issue with my initial composition model that results in an unphysical viscosity structure.

I am using ASCII data from the CRUST1.0 model to define the initial composition field. The original data has a resolution of 1° × 1° and my model size is 6° x 6° x 660 km, with a resolution of 0.25° x 0.25° x 10 km. The final composition results are shown in the first figure.

The red and blue represent the upper crust and lower crust, respectively. They have different rheological parameters and result in an unusual viscosity structure. Below is the viscosity structure picture.

This unusual viscosity structure can also lead to unphysical features in other fields. For example, the stress invariant shown in the third panel exhibits a rectangular anomaly, which I have highlighted with a red box.

My goal is to use real crustal layering structures, and I tried to improve the resolution of the model, but the results seemed to be the same.

By the way, the first and second images were taken at z = 640 km, and the third image was taken at z = 645 km. Here are some of my settings

  set Model name = visco plastic

  subsection Visco Plastic

    # Reference temperature and viscosity
    set Reference temperature = 273
    #set Reference viscosity = 1e22

    # 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-19
    set Reference strain rate = 1.e-16

    # Limit the viscosity with minimum and maximum values
    set Minimum viscosity = 1e19
    set Maximum viscosity = 1e24

    set Define thermal conductivities = true
    set Densities = background:3300, upper_crust:2600, lower_crust:2900
    set Heat capacities = 1300
    set Thermal conductivities = 3.6, 3, 3.2
    # density heteogenity for curst is layered but not affected by temperature because the density is already determined by CRUST1.0
    # density heteogenity for mantle is determined by inhomogeneous temperature field
    set Thermal expansivities = background:3e-5, upper_crust:0.0, lower_crust:0.0

    # 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 = composite

    set Grain size = 1e-3
    set Grain size exponents for diffusion creep =     3
    set Prefactors for diffusion creep            = background:2.37e-15, upper_crust:5.00e-51, lower_crust:5.00e-51
    set Activation energies for diffusion creep  =   background:3.75e+05, upper_crust:0.00e+00, lower_crust:0.00e+00
    set Activation volumes for diffusion creep   =   background:1.00e-05, upper_crust:0.00e+00, lower_crust:0.00e+00

    # 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)
    #    Gleason and Tullis (1995), Tectonophysics, v.472, p.213–225.
    #    " A flow law for dislocation creep of quartz aggregates determined with the molten salt cell"
    # Note that the viscous pre-factors below are scaled to plane strain from unixial strain experiments.
    set Prefactors for dislocation creep           =  background:6.52e-16, upper_crust:8.57e-28, lower_crust:7.13e-18
    set Stress exponents for dislocation creep    =  background:3.5, upper_crust:4.0, lower_crust:3.0
    set Activation energies for dislocation creep =  background:5.30e+05, upper_crust:2.23e+05, lower_crust:3.45e+05
    set Activation volumes for dislocation creep  =  background:1.80e-05, upper_crust:0.00e+00, lower_crust:0.00e+00

  end
end

Thank you in advance for your insights and suggestions!

Best regards,
Haolin

[bangerth]

@XHL-CUG This is really not my field, and so I can’t help with the actual question. But can I ask you to more clearly state what your question is? The way I read it is that you say your simulations “result in an unusual viscosity structure” (without actually saying what structure you are referring to, and why it is unusual). Underlying this question is, I think, that you have specific *expectations* what you think should happen, and you are assuming your reader to know what your expectations are and why. Be clear in your communications: “I expect that X would happen (for reasons A, B, C), but Y happens instead.”

Best
WB

[XHL-CUG]

@bangerth Thank you for your feedback.

My main question is: Why does Figure 2 show a rectangular region of high viscosity that is inconsistent with the surrounding viscosity field? What I would expect instead is a smooth viscosity field

To investigate this, I considered whether the issue might come from my initial compositional field (Figure 1), where the red region represents upper crust and the blue region represents lower crust. The CRUST1.0 input data I used has a resolution of 1° × 1° × 1 km, and the two compositional layers have different rheological parameters. This means that at certain depths, the model does contain a rectangular lower-crust region simply because the input data has that geometry.

After thinking more about it, I now suspect that this behavior may be unavoidable:
the rectangular high-viscosity feature in Figure 2 might simply reflect the rectangular lower-crust patch at that depth in Figure 1.

In other words, the sharp viscosity discontinuity may be a direct consequence of the blocky geometry of the CRUST1.0 compositional field. So this appears to be a data resolution issue rather than an issue with ASPECT.

Thank you again for your helpful feedback.

Best regards,
Haolin