Abnormal sediment_1 Distribution and Model Crash During Plate Subduction with LAB Thermal Adjustments

ASPECT
1

Hi everyone,

I’m encountering an unexpected issue with the sediment_1 component field in my ASPECT subduction models under Archean (Ar) and Phanerozoic (Ph) thermal regimes. When increasing the LAB temperature or reducing lithospheric thickness, half of the sediment develops a wedge-shaped morphology (see attached please), which destabilizes the crustal layer and ultimately causes the model to terminate with errors after significant slab tilting. Additionally, sporadic sediment_1 accumulation occurs intermittently on the opposite side of the domain. I attempted disabling the shared plugin, incrementally adjusting thermal parameters and lithospheric thickness, and deactivating FastScape, but these did not resolve the issue. Could anyone advise on the potential causes and solutions? Thank you for your expertise!

Best,
Noriel

Nomal one:

image
image780×432 9.09 KB

Abnormal one:
image
image828×464 37.6 KB

image
image1097×642 22.6 KB

It’s one of my input files:

Ar3-2.prm (20.4 KB)

###---------------------------------------------------------------------------Temperature------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Temperature boundary conditions
# Top and bottom (fixed) temperatures are consistent with the initial temperature field
# Note that while temperatures are specified for the model sides, these values are not used as the sides are not specified "Fixed temperature boundaries".  Rather, these boundaries are insulating (zero net heat flux).

subsection Boundary temperature model
  set Fixed temperature boundary indicators = bottom, top
  set List of model names = box

  subsection Box
    set Bottom temperature = 1825 
    set Top temperature    =  273
  end
end

subsection Initial temperature model
  set Model name = function

  subsection Function
    set Variable names = x,y
    set Function constants = h=600e3, ts1=273, ts2=723, ts3=993,\
                             A1=1.8e-6, A2=0.45e-6, A3=0.0, \
                             k1=2.5, k2=2.5, k3=2.5, \
                             qs1=0.074, qs2=0.038, qs3=0.020
	set Function expression = if((h-y) <= 20.e3, \
	    ts1 + (qs1/k1)*(h-y) - (A1*(h-y)*(h-y))/(2.0*k1), \
	    if((h-y) > 20.e3 && (h-y) <= 40.e3, \
	        ts2 + (qs2/k2)*(h-y-20.e3) - (A2*(h-y-20.e3)*(h-y-20.e3))/(2.0*k2), \
	        if((h-y) > 40.e3 && (h-y) <= 120.e3, \
	            ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3), \
	            if((h-y) > 120.e3, \
	                ts3 + (qs3/k3)*(120.e3-40.e3) - (A3*(120.e3-40.e3)*(120.e3-40.e3))/(2.0*k3) + 0.4*((h-y)-120.e3)/1.e3, 273 ))))
  end
end

# Constant internal heat production values (W/m^3) for background material
# and compositional fields.
subsection Heating model
  set List of model names = compositional heating

  subsection Compositional heating
    set Use compositional field for heat production averaging = 1, 0, 0, 1, 1, 1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0
    set Compositional heating values = 0.0, 0.0, 0.0, 1.8e-6, 0.45e-6, 0.0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0
  end
end
2

Hi @Noriel,

Thanks for posting the questions to the forum. I unfortunately don’t have any immediate insights into what exactly is happening.

Do you have any subduction models that are behaving in a manner that you expect?

It will be much easier to diagnose the issue if you can show how changing one or two parameters from a “well-behaved” model results in unphysical behavior.

Cheers,
John

3

Hi John，

I’m very sorry for the late reply—I took a few days off due to a fever and have now recovered. I wanted to show you the results of the “well-behaved” model and explain their differences.

First, in the “well-behaved” model, I set a lithospheric thickness of 140 km and a LAB temperature of 1633 K, which produced normal and expected results. Then, I adjusted the lithospheric thickness to 120 km and modified parameters like heat flow, heat generation, and LAB temperature to create a hotter thermal structure. However, at the start of the model run, sediments did not deposit at the surface as intended (I used Fastscape to control surface processes); instead, they appeared within the crust. This internal placement of sediments (rather than at the surface where they should be) caused uncontrollable errors, such as stratigraphic tilting, prolonged model runtime, and failure to form subduction.

In short, since I only modified the thermal structure and lithospheric thickness, I suspect other overlooked parameters (like the initial compositional field) might be responsible for this anomaly. Recently, I tried a “mid-ocean ridge” model by adjusting spreading velocity, successfully achieving a thinner lithosphere and higher LAB temperature, but I’m still curious about this error.

Thank you very much for your reply and insightful suggestions! :blush:

Best，
Noriel

Nomal one:

subsection Initial temperature model
  set Model name = function

  subsection Function
    set Variable names = x,y
    set Function constants = h=600e3, ts1=273, ts2=657, ts3=897,\
                             A1=1.2e-6, A2=0.30e-6, A3=0.0, \
                             k1=2.5, k2=2.5, k3=2.5, \
                             qs1=0.059, qs2=0.035, qs3=0.018
	set Function expression = if((h-y) <= 20.e3, \
	    ts1 + (qs1/k1)*(h-y) - (A1*(h-y)*(h-y))/(2.0*k1), \
	    if((h-y) > 20.e3 && (h-y) <= 40.e3, \
	        ts2 + (qs2/k2)*(h-y-20.e3) - (A2*(h-y-20.e3)*(h-y-20.e3))/(2.0*k2), \
	        if((h-y) > 40.e3 && (h-y) <= 140.e3, \
	            ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3), \
	            if((h-y) > 140.e3, \
	                ts3 + (qs3/k3)*(140.e3-40.e3) - (A3*(140.e3-40.e3)*(140.e3-40.e3))/(2.0*k3) + 0.4*((h-y)-140.e3)/1.e3, 273 ))))
  end
end

subsection Heating model
  set List of model names = compositional heating

  subsection Compositional heating
    set Use compositional field for heat production averaging = 1, 0, 0, 1, 1, 1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0
    set Compositional heating values = 0.0, 0.0, 0.0, 1.2e-6, 0.3e-6, 0.0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0
  end
end

Abnormal one:


subsection Initial temperature model
  set Model name = function

  subsection Function
    set Variable names = x,y
    set Function constants = h=600e3, ts1=273, ts2=723, ts3=993,\
                             A1=1.8e-6, A2=0.45e-6, A3=0.0, \
                             k1=2.5, k2=2.5, k3=2.5, \
                             qs1=0.074, qs2=0.038, qs3=0.020
	set Function expression = if((h-y) <= 20.e3, \
	    ts1 + (qs1/k1)*(h-y) - (A1*(h-y)*(h-y))/(2.0*k1), \
	    if((h-y) > 20.e3 && (h-y) <= 40.e3, \
	        ts2 + (qs2/k2)*(h-y-20.e3) - (A2*(h-y-20.e3)*(h-y-20.e3))/(2.0*k2), \
	        if((h-y) > 40.e3 && (h-y) <= 120.e3, \
	            ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3), \
	            if((h-y) > 120.e3, \
	                ts3 + (qs3/k3)*(120.e3-40.e3) - (A3*(120.e3-40.e3)*(120.e3-40.e3))/(2.0*k3) + 0.4*((h-y)-120.e3)/1.e3, 273 ))))
  end
end

subsection Heating model
  set List of model names = compositional heating

  subsection Compositional heating
    set Use compositional field for heat production averaging = 1, 0, 0, 1, 1, 1,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0
    set Compositional heating values = 0.0, 0.0, 0.0, 1.8e-6, 0.45e-6, 0.0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0
  end
end

The parameters I suspect：

subsection Boundary composition model
  set Model name = function
  set Fixed composition boundary indicators = top, bottom
  set Allow fixed composition on outflow boundaries = true
  
  subsection Function
    set Coordinate system   = cartesian
    set Variable names      = x,y,t
    set Function constants  = tsed=2e6 # duration of sedimentation switch time intervall 
    set Function expression           = 0; \
                                        0; \
                                        0; \
		           0; \
                                        0; \
                                        if(y>20e3&& t>0, 1, 0); \
                                        0; \
                                        if(y>20e3 && sign( sin( t/tsed * pi) ) == 1,1,0); \
                                        if(y>20e3 && sign( sin( t/tsed * pi) ) == -1,1,0); \
                                        0; \
                                        0; \
                                        0; \
                                        0; \
                                        0; \
                                        0; \
                                        0; \
                                        0  
  end
end
4

Hi @Noriel,

Thanks for posting the prm snippets and I agree the issue is likely related to the Initial composition and Boundary composition models.

Can you post both of those sections for both models?

Cheers,
John

5

Hi John,

This reply is my original input file. Please let me know if you need any further information :blush:

This is the failed model
Ar3-2.prm (20.4 KB)

And this is the successful model:
Ph-2.prm (19.9 KB)

Best，
Noriel

6

Hi @Noriel,

I think the issue does indeed lie in the Boundary composition model.

I suggest the following changes:

  1. Refer to interface locations in terms of height rather than depth (i.e., y>20e3 → y>580e3)
  2. Set the function expression equal to 0 for all fields not associated with sediments.
  3. Only have one compositional field representing sediment.

I’ve included an example PRM snippet below from a model my group ran recently with fastscape coupling, which will hopefully give context for my suggestions above.

Cheers,
John

subsection Compositional fields
  set Number of fields = 10
  set Names of fields  = plastic_strain, \
                         noninitial_plastic_strain, \
                         viscous_strain, \
                         crust_upper, \
                         crust_lower, \
                         mantle_lithosphere, \
                         asthenosphere, \
                         sediment, \
                         sediment_age, \
                         deposition_depth 
  set Types of fields = strain, \
                        strain, \
                        strain, \
                        chemical composition, \
                        chemical composition, \
                        chemical composition, \
                        chemical composition, \
                        chemical composition, \
                        generic, \
                        generic
  set Compositional field methods = particles
  set Mapped particle properties  = plastic_strain: plastic_strain, \
                                    noninitial_plastic_strain: noninitial_plastic_strain, \
                                    viscous_strain: viscous_strain, \
                                    crust_upper: initial crust_upper, \
                                    crust_lower: initial crust_lower, \
                                    mantle_lithosphere: initial mantle_lithosphere, \
                                    asthenosphere: initial asthenosphere, \
                                    sediment: initial sediment, \
                                    sediment_age: initial sediment_age, \
                                    deposition_depth: initial deposition_depth

end

# Initial values of different compositional fields
# The upper crust (20 km thick), lower crust (20 km thick)
# and mantle (60 km thick) are continuous horizontal layers
# of constant thickness. The viscoelastic stresses and 
# non initial plastic strain are set to 0, while the
# initial plastic strain is randomized between 0.5 and 1.5.
subsection Initial composition model
  set List of model names = function, world builder
  set List of model operators = add, add

  subsection Function
    set Variable names      = x,y
    set Function expression = if(x>140.e3 && x<260.e3 && y>340.e3, 0.5 + rand_seed(0.1), 0); \
                              0; \
                              0; \
                              if(y>=380.e3, 1, 0); \
                              if(y<380.e3 && y>=360.e3, 1, 0); \
                              if(y<360.e3 && y>=320.e3, 1, 0); \
                              if(y<320.e3 && y>=0.e3, 1, 0); \
                              0; \
                              0; \
                              0;

  end
end

subsection Boundary composition model
  set Model name =  function
  set Fixed composition boundary indicators = top, bottom
  set Allow fixed composition on outflow boundaries = true

  subsection Function
    set Coordinate system   = cartesian
    set Variable names      = x,y,t
    set Function constants  = DZ = 400e3, SeaLevel = -100  #  Attention: 2D setup only. Make sure to update these parameters according to box size and fastscape parameters. 
    set Function expression = \
                             0; \
                             0; \
                             0; \
                             0; \
                             0; \
                             0; \
                             0; \
                             if(y>390e3, t/1e6, 0); \
                             if(y>390e3, DZ - SeaLevel -y, 0); \
                             if(y>390e3, 1, 0); 
  end
end