Element shape for "Box with lithosphere boundary indicators" geometry model

ASPECT
1

Good afternoon,

I’m running 2D subduction models in ASPECT, where I’m using the “box with lithosphere boundary indicators” geometry model to impose horizontal velocities at the side of the model domain and at lithospheric depths (<= 100 km). Because I have an elongated domain (aspect ratio 1:6; x,y = 5400 km, 900 km), and am using the lithospheric indicators, my problem is that I can’t seem to get square finite elements throughout the entire domain without introducing issues into the subduction model:

If I generate the mesh by choosing 1 repetition along the Y axis, I get square elements in all the domain except in the lithosphere, where the elements are rectangular and with aspect ratio 1:6. In particular:

subsection Geometry model
set Model name = box with lithosphere boundary indicators

subsection Box with lithosphere boundary indicators
set X extent = 5400.e3
set Y extent = 900.e3
set X repetitions = 6
set Y repetitions = 1
set Lithospheric thickness = 100.e3
set Y repetitions lithosphere = 1
end
end

If I use the lithospheric thickness as a scaling factor, and have 9 repetitions along y and 54 along x, i.e. replacing the corresponding two lines with the following:

set X repetitions = 54
set Y repetitions = 9

I do get the square elements, but inflow boundary conditions do not work as expected (horizontal velocity is no longer imposed at just lithospheric depths) . The prm files I used for these two tests (attached) are identical except for the X, Y repetitions I reported above.
rectangulat_elements.prm (9.8 KB)
square_elements.prm (9.8 KB)

Do you know what the issue is here? Any help would be highly appreciated.

Thanks,
Valeria Turino

2

Hi Valeria,

Apologies, I forgot to reply to your message when I first saw it. I have not had a chance to rerun your models, but below is an initial response.

The presence of distortion in your first setup makes sense as you should end up with one vertical element in the lower domain and one vertical element in the lithosphere. However, in the lithosphere I would have predicted a 9x1 element ratio.

For reference, In 2D models the parameter Y repetitions refers to the number of vertical elements in the lower domain, and the parameter Y repetitions lithosphere refers to the number of vertical elements in the lithosphere (see documentation here)

For the second scenario, I think achieving exactly square elements would require the following parameters:

subsection Box with lithosphere boundary indicators
  set X extent = 5400.e3
  set Y extent = 900.e3
  set X repetitions = 54
  set Y repetitions = 8
  set Lithospheric thickness = 100.e3
  set Y repetitions lithosphere = 1
end

The key change above is set Y repetitions = 8, as the whole domain is 900 km in the y-direction, but you only have 800 km in the vertical direction in the lower region. As noted further down in the parameter documentation, if a cell edge does not lie exactly on the lithosphere boundary, the applied velocities will not be exactly at that depth.

I wonder if the latter is causing the reported issues?:

I do get the square elements, but inflow boundary conditions do not work as expected (horizontal velocity is no longer imposed at just lithospheric depths)

Can you retry with the proposed adjustment above and see if that solves the issue?

Cheers,
John

3

Hi John,
Thank you for your reply!

I tried re-running the model following your suggestion, but I still get the same issue I was having previously with the velocity field.

I’m attaching the new parameter file with the modified Y repetitions as suggested.

Thanks,
Valeria
square_elements_new.prm (9.8 KB)

4

Hi Valeria,

Thanks for trying out that suggestion! Can you send over images detailing the mismatch between the expected and observed velocity field, as well as the mesh?

I will not be at a location where I can run the model for a bit, but will be able to do so tomorrow.

Cheers,
John

5

Hi John and Valeria,

I can do that, since I was also trying to figure out what the problem is, but I still don’t have a good idea.

Screenshot from 2023-06-16 12-48-56
Screenshot from 2023-06-16 12-48-561844×1116 346 KB

This picture shows the velocity in X direction (colorscale between 4 and 6 cm/yr), and the X velocity is supposed to be 5 cm/yr in the lithosphere. You can see that the value is actually correct in the lithosphere, and it’s different from 5 cm/yr below that (so it’s not like the wrong value is being prescribed below the lithosphere), but I do not understand why the two models are so different. The resolution is the only thing that should be different between the two.

Juliane

6

Hi Juliane and Valeria,

Thanks for running those models and posting the images! To confirm, the top image is from square_elements.prm and the second images is from square_elements_new.prm?

I’m not sure either why the two images are so different. Not having the elements align on the lithosphere boundary must be making a large difference?

I’m also confused by the values of the x velocities.

The boundary velocity section contains:

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

  subsection Function
    set Function constants  = cm=100.0 #*3600.0*24.0*365.25
    set Function expression = if(x==0, 0.05, \
                                 if(x==5400.e3, -0.04, 0));\
                              0
    set Variable names      = x,z
  end
end

From the above equation, shouldn’t the x velocity on the right boundary should be negative 0.04?

Valeria - can you try rerunning the simulation with the following boundary velocity settings?

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

subsection Function
set Function constants = cm=100.0 #3600.024.0*365.25
set Function expression = if(x<2700.e3, 0.05, -0.05); 0
set Variable names = x,z
end
end

The above conditions should give a symmetric velocity field, and if that works then we can start re-adding complexity piece by piece.

Cheers,
John

7

Hi John,
No, I was running the

and I was just looking at the velocity on the left boundary, since that seemed to be the place where the difference occurs (below the lithosphere, the boundary condition is a traction boundary condition for prescribed lithostatic pressure, so if something is different in the model, it would make sense that you would get a difference there). I am just not sure what causes the difference.

Juliane

8

Hi John and Juliane,
thanks for all the help!

I tried running a model with John’s latest suggestion and I get a better velocity field now, but not the same as the original rectangular elements model.
Using the same scale Juliane had:

image
image1918×586 10.9 KB

Using a symmetric scale:

image
image1915×616 35.4 KB

Also, the model runs until timestep 117288, then the iterative advection solver stops converging. I’m guessing it’s a conservation of mass issue because none of the boundaries are open with the new BCs?

Here’s the new parameter file:
square_elements_symmBCs.prm (9.7 KB)

Thanks again!
Valeria