How to set a continental plate with different depth?

Dear friend: I installed the WB successfully, thank you. I'm trying to build a subduction model with a thin oceanic plate and a thick continental plate. The subducting slab with a small dip is too thin and it through the continental plate in the middle. I tried to build several continental plates with different depth to avoid this but the aspect result is not convergent because of the existing null space.(I built several mantle layer blocks to insure no null space and the result is still error.). Whether it allowes a continental plate with different depth or not? Or how to avoid the thin slab through the thick continental plate?

{ "version":"0.3", "cross section":[[0,0],[800e3,0]], "features": [

    {
    "model":"continental plate","name":"South America1","max depth":100e3,
        "coordinates":[[550e3,-1],[800e3,-1],[800e3,1],[550e3,1]],
        "temperature models":[{"model":"linear","max depth":100e3,"bottom temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[1],"max depth":20e3},
                    {"model":"uniform","compositions":[2],"min depth":20e3,"max depth":40e3},
                    {"model":"uniform","compositions":[3],"min depth":40e3,"max depth":100e3}]
    },

    {
    "model":"continental plate","name":"South America2","max depth":90e3,
        "coordinates":[[500e3,-1],[550e3,-1],[550e3,1],[500e3,1]],
        "temperature models":[{"model":"linear","max depth":90e3,"bottom temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[1],"max depth":20e3},
                    {"model":"uniform","compositions":[2],"min depth":20e3,"max depth":40e3},
                    {"model":"uniform","compositions":[3],"min depth":40e3,"max depth":90e3}]
    },

    {
    "model":"continental plate","name":"South America3","max depth":60e3,
        "coordinates":[[410e3,-1],[500e3,-1],[500e3,1],[410e3,1]],
        "temperature models":[{"model":"linear","max depth":60e3,"bottom temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[1],"max depth":20e3},
                    {"model":"uniform","compositions":[2],"min depth":20e3,"max depth":40e3},
                    {"model":"uniform","compositions":[3],"min depth":40e3,"max depth":60e3}]
    },

    {
    "model":"continental plate","name":"Weak Zone","max depth":60e3,
        "coordinates":[[400e3,-1],[410e3,-1],[410e3,1],[400e3,1]],
        "temperature models":[{"model":"linear","max depth":60e3,"bottom temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[0],"max depth":60e3}]
    },

    {
    "model":"oceanic plate","name":"Nazca","min depth":0e3,"max depth":55e3,"coordinates":[[0,-1],[400e3,-1],[400e3,1],[0,1]],
        "temperature models":[{"model":"plate model","min depth":0,"max depth":55e3,"spreading velocity":0.01,"ridge coordinates":[[200e3,-1],[200e3,1]],"top temperature":273,"bottom temperature":1573}],
        "composition models":[
                    {"model":"uniform","compositions":[4],"max depth":4e3},
                    {"model":"uniform","compositions":[5],"min depth":4e3,"max depth":12e3},
                    {"model":"uniform","compositions":[6],"min depth":12e3,"max depth":55e3}
                     ]
    },

    {
    "model":"mantle layer","name":"Asthonisphere1","min depth":55e3,"max depth":400e3,
        "coordinates":[[0,-1],[400e3,-1],[400e3,1],[0,1]],
        "temperature models":[{"model":"linear","min depth":55e3,"max depth":400e3,"top temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[7]}]
    },

    {
    "model":"mantle layer","name":"Asthonisphere2","min depth":60e3,"max depth":400e3,
        "coordinates":[[400e3,-1],[500e3,-1],[500e3,1],[400e3,1]],
        "temperature models":[{"model":"linear","min depth":60e3,"max depth":400e3,"top temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[7]}]
    },

    {
    "model":"mantle layer","name":"Asthonisphere3","min depth":90e3,"max depth":400e3,
        "coordinates":[[500e3,-1],[550e3,-1],[550e3,1],[500e3,1]],
        "temperature models":[{"model":"linear","min depth":90e3,"max depth":400e3,"top temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[7]}]
    },

    {
    "model":"mantle layer","name":"Asthonisphere4","min depth":100e3,"max depth":400e3,
        "coordinates":[[550e3,-1],[800e3,-1],[800e3,1],[550e3,1]],
        "temperature models":[{"model":"linear","min depth":100e3,"max depth":400e3,"top temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[7]}]
    },

    {
    "model":"subducting plate","name":"Nazca Slab","coordinates":[[400e3,-1],[400e3,1]],"dip point":[750e3,400e3],
        "segments":[{"length":50e3,"thickness":[55e3],"angle":[0,15]},
        {"length":200e3,"thickness":[55e3],"angle":[15,40]}],
        "temperature models":[{"model":"plate model","plate velocity":0.078}],
        "composition models":[{"model":"uniform","compositions":[4],"max distance slab top":4e3},
        {"model":"uniform","compositions":[5],"min distance slab top":4e3,"max distance slab top":12e3},
        {"model":"uniform","compositions":[6],"min distance slab top":12e3,"max distance slab top":55e3}]
    }

]}

Am I understanding you correctly that you are trying to run this through the geodyanmics code ASPECT, and that ASPECT crashes with a convergence error?

Or is there something wrong with the initial conditions which are being generated? If so, can you attach a picture of what is generated and explain why it is not what you expect/want?

test3 test2

Yes. That's what you mean. Aspect run the upper model successfully and failed in running the lower one.

Thanks for the pictures. Without the aspect input file or a detailed error output it is hard to help you with debugging the problem, but I have a few comments about the world builder fil e and the output:

Why is the South America3 continental plate 600 km thick? Is that intentional or a typo?
A way to have a mantle temperature below the slab is to make the slab thicker and divide it up in an upper part containing the plate model temperature model and a lower part containing an adibatic temperterature model (like you do with the compositional model).
oceanic plate thickness seem so be different (thicker) than the thickness of the subducting plate. You might want to look into whether this is what you want in your case.

If that doesn't help, can you provide a bit more information on the model and/or error output?

Thank you. The 60km thick continental plate(South America3) was intentional to have a mantle temperature below the slab. I change the thickness of the slab to make it as thick as the oceanic plate in a new large size model, but the mistake is still exist. The mistake may not exist when building the model as the upper picture. I paste my aspect input file and the world builder file below and also the error. And How can I set the WB to make the temperature field continuous in the yellow circle？

T20ma aspect input file: `set Output directory = /home/ubuntu/changs/data set Use years in output instead of seconds=true set World builder file = /home/ubuntu/aspect/input/5ridgesub-wb.wb set Dimension = 2 set CFL number = 0.5 set End time = 3e6 set Pressure normalization =no set Maximum time step = 1e10 set Max nonlinear iterations = 100000000 set Nonlinear solver tolerance = 1e-5

subsection Solver parameters subsection Stokes solver parameters set Maximum number of expensive Stokes solver steps = 10000000 end end

subsection Geometry model set Model name = box subsection Box set X extent = 2000e3 set Y extent = 660e3 set X repetitions = 6 set Y repetitions = 2 end end

subsection Material model set Model name = visco plastic subsection Visco Plastic set Thermal diffusivities=9.89e-7,1.21e-6,1.21e-6,1.15e-6,9.87e-7,1.21e-6,1.17e-6,1.17e-6 set Specific heats = 750 set Densities = 3300,3300,2800,2900,3300,3300,3300,3300 set Thermal expansivities =2e-5 set Angles of internal friction=20,2,20,20,20,5,10,20 set Cohesions=20e6,1e6,20e6,20e6,20e6,1e6,1e6,20e6 set Prefactors for diffusion creep=2.37e-15,1e-50,1e-50,1e-50,2.37e-15,1e-50,1e-50,2.37e-15 set Prefactors for dislocation creep=6.52e-15,1.12e-10,8.57e-28,7.13e-18,6.52e-16,1.12e-10,1.12e-10,6.52e-16 set Stress exponents for dislocation creep=3.5,3.4,4.0,3,3.5,3.4,3.4,3.5 set Activation energies for diffusion creep=375e3,0,0,0,375e3,0,0,375e3 set Activation energies for dislocation creep=530e3,497e3,223e3,345e3,530e3,497e3,497e3,530e3 set Activation volumes for diffusion creep=4e-6,0,0,0,4e-6,0,0,4e-6 set Activation volumes for dislocation creep=13e-6,0,0,0,13e-6,0,0,13e-6 set Grain size exponents for diffusion creep=3,1,1,1,3,1,1,3 end end

subsection Compositional fields set Number of fields = 7 end

subsection Initial composition model set Model name = world builder end

subsection Gravity model set Model name = vertical subsection Vertical set Magnitude = 9.8 end end

subsection Mesh refinement set Additional refinement times =0 set Initial adaptive refinement =2 set Initial global refinement =3 set Refinement fraction = 0.7 set Coarsening fraction = 0.05 set Time steps between mesh refinement =10 end

subsection Initial temperature model set Model name = world builder end

subsection Boundary temperature model set Fixed temperature boundary indicators=2,3 set Model name = box subsection Box set Bottom temperature = 1713 set Top temperature = 273 end end

subsection Boundary velocity model set Prescribed velocity boundary indicators = left:function set Tangential velocity boundary indicators = right,bottom end

subsection Mesh deformation set Mesh deformation boundary indicators=top:free surface subsection Free surface set Free surface stabilization theta=0.5 end end

subsection Boundary velocity model subsection Function set Variable names = x,y set Function constants = Va=0.078 set Function expression = if(y>560e3,Va,0);0; end end

subsection Postprocess set List of postprocessors = visualization,composition statistics, temperature statistics,topography subsection Visualization set Time between graphical output = 1e5 set List of output variables = density, temperature anomaly set Interpolate output = true end end World builder file { "version":"0.3", "cross section":[[0,0],[2000e3,0]], "features": [

    {
    "model":"mantle layer","name":"Asthonisphere","min depth":63e3,"max depth":660e3,
        "coordinates":[[0,-1],[1150e3,-1],[1150e3,1],[0,1]],
        "temperature models":[{"model":"linear","min depth":63e3,"max depth":660e3,"top temperature":1563}],
        "composition models":[{"model":"uniform","compositions":[7]}]
    },

    {
    "model":"mantle layer","name":"Asthonisphere","min depth":100e3,"max depth":660e3,
        "coordinates":[[1150e3,-1],[2000e3,-1],[2000e3,1],[1150e3,1]],
        "temperature models":[{"model":"linear","min depth":100e3,"max depth":660e3,"top temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[7]}]
    },

    {
    "model":"continental plate","name":"South America","max depth":100e3,
        "coordinates":[[1150e3,-1],[2000e3,-1],[2000e3,1],[1150e3,1]],
        "temperature models":[{"model":"linear","max depth":100e3,"bottom temperature":1573}],
        "composition models":[{"model":"uniform","compositions":[1],"max depth":20e3},
                    {"model":"uniform","compositions":[2],"min depth":20e3,"max depth":40e3},
                    {"model":"uniform","compositions":[3],"min depth":40e3,"max depth":100e3}]
    },

    {
    "model":"continental plate","name":"Weak Zone","max depth":63e3,
        "coordinates":[[1000e3,-1],[1150e3,-1],[1150e3,1],[1000e3,1]],
        "temperature models":[{"model":"linear","max depth":63e3,"bottom temperature":1563}],
        "composition models":[{"model":"uniform","compositions":[0],"max depth":100e3}]
    },

    {
    "model":"oceanic plate","name":"Nazca","min depth":0e3,"max depth":63e3,"coordinates":[[0,-1],[1000e3,-1],[1000e3,1],[0,1]],
        "temperature models":[{"model":"plate model","min depth":0,"max depth":63e3,"spreading velocity":0.031,"ridge coordinates":[[380e3,-1],[380e3,1]],"top temperature":273,"bottom temperature":1563}],
        "composition models":[
                    {"model":"uniform","compositions":[4],"max depth":4e3},
                    {"model":"uniform","compositions":[5],"min depth":4e3,"max depth":12e3},
                    {"model":"uniform","compositions":[6],"min depth":12e3,"max depth":63e3}
                     ]
    },

    {
    "model":"subducting plate","name":"Nazca Slab","coordinates":[[1000e3,-1],[1000e3,1]],"dip point":[1800e3,660e3],
        "segments":[{"length":120e3,"thickness":[63e3,68e3],"angle":[0,10]},
        {"length":190e3,"thickness":[68e3,77e3],"angle":[10,30]},
        {"length":100e3,"thickness":[77e3],"angle":[30,40]}],
        "temperature models":[{"model":"plate model","plate velocity":0.078}],
        "composition models":[{"model":"uniform","compositions":[4],"max distance slab top":4e3},
        {"model":"uniform","compositions":[5],"min distance slab top":4e3,"max distance slab top":12e3},
        {"model":"uniform","compositions":[6],"min distance slab top":12e3,"max distance slab top":77e3}]
`

error error

The 60km thick continental plate(South America3) was intentional to have a mantle temperature below the slab.

There is a background temperature by default, which is an adiabatic temperature. So there is no need to do this. If you want to have a different temperature for the mantle, I would strongly suggest using the mantle layer feature.

I change the thickness of the slab to make it as thick as the oceanic plate in a new large size model, but the mistake is still exist. The mistake may not exist when building the model as the upper picture.

It does look a bit better now, but there is a fundamental problem you are encountering. The subducting plate temperature model plate model is based on the McKenzie 1970 paper (I would like to rename this temperature model to McKenzie 1970 to make it more clear). It assumes that the oceanic plate is old and therefore has a linear temperature profile. So if you want a perfect connection you either have to

put the ridge of the oceanic plate further away so that the oceanic plate temperature has become linear by the time it reaches the trench.
come up with a temperature model which doesn't assume the incoming oceanic plate to be very old, or wait for this the be implemented into the world builder (I know that there is more interest in this, and that others are planning to work on this). If you have a temperature model you would like to implement for the subducting plate I can help you with that.

But if the disconnect is not too bad, it may not influence the part of the model you are interested in, but you would need to verify that.

error

This is more of an aspect question than a world builder question since it only happens at the 10th timestep, but I might be able to give you a few tips to help you debug the problem. I ran the problem myself and noticed that it only crashed when the timestep fluctuating a lot (0-1: 29, 1-2: 18, 2-3: 27, 3-4: 40, 4-5: 42, 5-6: 53, 6-7: 43, 7-8: 72, 8-9: 44, 9-10:196 yr). This could be an indication that you model is unstable (I am especially worried about the free surface in this case). My hypothesis would be that the free surface has a hard time to stabilize because there is a temperature difference going all the way through the mantle.

I would suggest you first run the model and output visualization output at every timestep. Look at the different fields for anomalies (is weird material/temperatures flowing into the boundaries) and look at how the velocity field changes between every timestep. Secondly try to run the model with the aspect parameter set Maximum relative increase in time step and set it to something like 50. I tried that out and it seems to be able to compute past 800 years without crashing.

Thank you very much, my problem is solved, I changed the CFL number to 0.1 and set the average model to harmonic average.

“The subducting plate temperature model plate model is based on the McKenzie 1970 paper; If you have a temperature model you would like to implement for the subducting plate I can help you with that.‘’

I think a half-space cooling model is not enough for many conditions, a model such as GDH1(stein and stein, 1992) or a plate cooling model mentioned in(F. Richards et al, 2018) is better.

I have a question. In WB, a subducting plate needs to set the plate velocity. So has the initial temperature of the slab taken into account the effect of mantle heating during its subduction? If it's been considered I‘d like to know how it's calculated. By the way, it would be even better if expressions of initial temperature fields could be explained in the manual.

Thank you. You've helped a lot.

Hey,

Thank you very much, my problem is solved, I changed the CFL number to 0.1 and set the average model to harmonic average.

It's great to hear that that problems is solved.

I think a half-space cooling model is not enough for many conditions, a model such as GDH1(stein and stein, 1992) or a plate cooling model mentioned in(F. Richards et al, 2018) is better.

Are you talking here about the oceanic plate or subducting plate temperature model? For either it is very easy to add new models. If you can write down the equations in terms of distance from the ridge and distance from the top of the plate for an oceanic plate and distance from the trench and distance from the top of the slab for the subducting plate, I can code them up for you, or, if you want, tell you how you can make the modules yourself and contribute them to the world builder :)

I have a question. In WB, a subducting plate needs to set the plate velocity. So has the initial temperature of the slab taken into account the effect of mantle heating during its subduction? If it's been considered I‘d like to know how it's calculated. By the way, it would be even better if expressions of initial temperature fields could be explained in the manual.

That is a good questions, which could indeed use some extra information in the manual. In short, the equation from the McKenzie 1970 paper assumes a linear temperature gradient at the trench and than heats that up as a halfspace. The only thing why the velocity it accounts for that the slab going down fast in the mantle has less time to heat up and will remain colder than a slab which subducts really slowly.

Thank you. You've helped a lot.

I am glad it has been helpful!

Hi,

For either it is very easy to add new models. If you can write down the equations in terms of distance from the ridge and distance from the top of the plate for an oceanic plate and distance from the trench and distance from the top of the slab for the subducting plate, I can code them up for you, or, if you want, tell you how you can make the modules yourself and contribute them to the world builder :)

Thank you. I paste the Plate Model(stein and stein, 1992) here. wb2 wb3

If it's convenient I'd like to know how to do it after the new model is established. ^_^

Another question here. The thickness of an oceanic plate changes with plate age(Picture below). I notice WB change the temperature field but didn't change the thickness of the ocean plate after setting up the ridge. Is it possible to change the plate thickness and the temperature field at the same time when set up the ridge coordinates and the spreading rate? For example, the plate boundary(thickness) is the Tm=1300 ℃ isothermal.

Thank you. I paste the Plate Model(stein and stein, 1992) here.

I think this is exactly what is in the oceanic plate feature temperature model called plate model, so it's already implemented :)

see https://github.com/GeodynamicWorldBuilder/WorldBuilder/blob/master/source/features/oceanic_plate_models/temperature/plate_model.cc#L172 for the implementation. If you need something else, let me know.

Another question here. The thickness of an oceanic plate changes with plate age(Picture below). I notice WB change the temperature field but didn't change the thickness of the ocean plate after setting up the ridge. Is it possible to change the plate thickness and the temperature field at the same time when set up the ridge coordinates and the spreading rate? For example, the plate boundary(thickness) is the Tm=1300 ℃ isothermal.

A more generalized version of this problem is quite high on my todo list, but I am not sure when I will be able to get around to implementing it. It is very much related to #132 in a way, because I need to find a good way to represent complicated 3D surfaces in the world builder internally. I would like to ask you the question of why you think you need that? Depending on your problems there may be good ways to solve your problem in a different way without changing code, or only having very minor changes.

If non of the tricks work for you and it is important to you and you need it quickly, then we can discuss it and I can help you implement it.

Thank you！ I just want to make the model more realistic ，not in a hurry.

You can now make layers change thickness spatially (see the new online manual for how to do this). Feel free to reopen if this doesn't answer your question.

GeodynamicWorldBuilder / WorldBuilder

How to set a continental plate with different depth? #207