Does ASPECT consider seawater density in calculating dynamic topography

Hi all,
I have been testing 3D cases with and without a free surface to see whether they yield the same dynamic topography. The difference between the two models are shown below:


The topography and dynamic topography (after removing the seafloor topography in a blank model without a plume) on the xz cross-section of the center of the plume at 200 Myr are shown below:

I wonder why there is a huge difference in topography and dynamic topography between the two cases. Also, may I ask: does ASPECT consider the seawater density above the free surface in the case with a free surface or does ASPECT consider the seawater density when calculating the dynamic topography in the model without a free surface? I would greatly appreciate it if anyone could help me answer these questions. Thank you very much!

Cheers,
Ziqi

Hi Ziqi,

I think we need a bit more information to clear up all your questions (ideally a parameter file for both models and screenshots of the solution). I would also recommend replacing the boundary numbers with boundary names so that it is easier to understand what happens at which boundary (the different behavior between free surface and no free surface at the right end of your plots – x=3000km – suggests something slightly weird is happening at that boundary are you prescribing a velocity or traction at that boundary?). Finally, can you clarify what you mean by topography and dynamic topography? I dont think you mean the output of the topography postprocessor, because that should be 0 and constant in the case of no free surface, and also not the output of the dynamic topography postprocessor, because that should be close to 0 in the case of a free surface

Also, may I ask: does ASPECT consider the seawater density above the free surface in the case with a free surface or does ASPECT consider the seawater density when calculating the dynamic topography in the model without a free surface?

This I can answer. The free surface assumption is that stress at the surface is exactly 0, i.e. no seewater above the surface is considered. You may be able to prescribe sea water if you combine the free surface with a traction boundary that prescribes a certain amount of pressure in dependence of depth, but I do not know if someone has tried and tested that case. The dynamic topography postprocessor has a parameter called Density above, which lets you choose which density needs to be replaced by the model material. It defaults to 0 (essentially air), but can be set to other values (e.g. 1000 for water).

Hope that helps a bit,
Rene

Hi Rene,
Thank you very much for your helpful reply!
Sorry, I wasn’t clear about my question earlier. So I tested two cases.
The first case has a free surface and the “box with lithosphere” geometry and boundary conditions are shown below:

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

  subsection Box with lithosphere boundary indicators
    set Box origin X coordinate = 0
    set Box origin Y coordinate = 0
    set X extent                = 3000e3
    set Y extent                = 500e3
    set Z extent                = 1000e3
    set X repetitions           = 60
    set Y repetitions           = 10
    set Z repetitions           = 20
    set Z repetitions lithosphere       = 1
    set Lithospheric thickness          = 50e3
  end
end
subsection Mesh deformation
  set Mesh deformation boundary indicators = top: free surface & diffusion
  set Additional tangential mesh velocity boundary indicators = left, right, front, back, lithosphere front, lithosphere back
  subsection Free surface
    set Surface velocity projection             = vertical
    set Free surface stabilization theta        = 1
  end
  subsection Diffusion
    set Hillslope transport coefficient  = 5e-8
  end
end
subsection Boundary traction model
  set Prescribed traction boundary indicators = left: initial lithostatic pressure, right: initial lithostatic pressure, front: initial lithostatic pressure, back: initial lithostatic pressure

  subsection Initial lithostatic pressure
    # We calculate the pressure profile at the right model boundary.
    set Representative point         = 3000e3, 500e3, 1000e3
  end
end

In this case, I implement the plume by adding an initial thermal anomaly and a prescribed velocity at the bottom boundary.
In the postprocessor, I output the topography shown below:


Then I ran a completely same case but without the plume (no initial thermal anomaly nor prescribed velocity at the bottom) and acquired the topography from the blank case.

I got the dynamic topography by subtracting the topography of the blank case from the first case.

My question for this case is: is the density above the free surface air or seawater (I add traction to the left, right, front and back boundaries)?
In the second case, I use the “box” geometry model and the same boundary conditions.

subsection Geometry model
  set Model name                = box

  subsection Box
    set Box origin X coordinate = 0
    set Box origin Y coordinate = 0
    set X extent                = 3000e3
    set Y extent                = 500e3
    set Z extent                = 1000e3
    set X repetitions           = 60
    set Y repetitions           = 10
    set Z repetitions           = 20
  end
end
subsection Boundary traction model
  set Prescribed traction boundary indicators = left: initial lithostatic pressure, right: initial lithostatic pressure, front: initial lithostatic pressure, back: initial lithostatic pressure
  subsection Initial lithostatic pressure
    # We calculate the pressure profile at the right model boundary.
    set Representative point         = 3000e3, 500e3, 1000e3
  end
end

In the postprocessor, I output the dynamic topography shown below:


Then I ran a completely same case but without the plume (no initial thermal anomaly nor prescribed velocity at the bottom) and acquired the dynamic topography from the blank case:

Finally, I got the real dynamic topography by subtracting the dynamic topography of the blank case:

For this case, I know the density above is 0 because that’s how I set it up in the dynamic topography postprocessor.
In the end, my ultimate question is which, of the above two dynamic topography, is the more accurate/robust dynamic topography so that I can use to compare with/match with the residual topography observed in nature?

Cheers,
Ziqi

Hi Ziqi,
Thanks for the additional details. I am still not sure I understand all the necessary details of your models, but let’s start with what I can answer:

For your case 1: The “density” above the free surface is not prescribed, but the effect of the boundary condition is similar to air (or more precisely to have no stress at the top boundary, so the most similar “density” would be vacuum). “Free surface” usually implies two things in geodynamical models: That the stress at the boundary is zero (which is a boundary condition for the Stokes equation) and that the boundary nodes of the mesh move freely following the velocity computed at the top boundary. This is what is implemented in ASPECT.

My questions for your case 1: I think something is wrong (or at least unusual) with your boundary conditions. You prescribed all side boundaries to be “Additional tangential mesh velocity boundary indicators” except for the boundaries lithosphere left and lithosphere right. That means the mesh at those boundaries cannot deform. These boundaries are what forces the topography to be zero at x=0 and x=3000km. Your topography could look very different if you added those boundaries to be tangential mesh velocity (but maybe you are aware of this already). Generally it looks like your model is loosing some volume, which you can see in the largely negative topography (for whichever reason, e.g. possibly cooling of the model, the vertical projection of surface velocity and the surface diffusion could also play a role). This motivates my next question:
In the subsection Boundary traction model you prescribe initial lithostatic pressure at all side boundaries except for the lithosphere (lithosphere left, lithosphere right, lithosphere front, lithosphere back). My question would be what boundary condition do you apply there? If you prescribe no boundary condition (neither traction nor velocity) that means material would pour without resistance out of the model at those boundaries (like an overflowing bucket), which could explain your negative topography.

For your case 2: This input file looks more like I expected, however I do not fully understand the output either. What is the cause of the strong negative topography at x=3000km? According to your boundary condition material should be able to enter/leave the domain at lithostatic conditions. Do you prescribe any velocities in your model except for the plume influx? Your model looks a bit like an aging plate, but how do you prescribe that motion? One idea I had was that the initial lithosphere boundary traction condition compute the traction profile at the initial time. If you run your model for 200 Myr the pressure at the boundary might change and lead to a flux out of the model. How does the negative topography at x=3000km evolve over time? And do you observe a significant average pressure change at x=3000km over time?

In general for identical setups (one with a free surface, one with the dynamic topography postprocessor) we would expect the two methods to deliver results that are qualitatively the same, with some deviations (<20-30%) in topography amplitude. The fact that your two models look so different points to some difference between the setups that you need to track down. It is not a fundamental difference in the ways to compute the topography at least. You could try to simplify your setups and make them more similar by using the same geometry model.

I hope that helps, for more advice I would need the full input files.

Best,
Rene

Hi Rene,
Thank you very much for your detailed and helpful reply.
First, to answer your questions: for both models, I implement a 2-cm-per-year plate velocity at the top, which is why I added the ‘box with lithosphere’ model in the first case because I can’t prescribe the plate velocity at the top boundary if it is free surface.
For the first case, thanks for pointing out why topography is forced to be zero at x=0 and x=3000km. The reason why I didn’t make lithosphere left and lithosphere right to ‘Additional tangential mesh velocity boundary’ is I implemented the prescribed velocity on those two boundaries.

subsection Boundary velocity model
set Prescribed velocity boundary indicators = bottom: function, lithosphere left: function, lithosphere right : function, front y: function, back y: function, lithosphere front y: function, lithosphere back y: function
  subsection Function
    set Variable names          = x,y,z
    set Function constants      = vx_top=0.02, vx_bot=0.0, vy_top=0, vy_bot=0, vz_bot=0.0, vz_top=0.0, x0=1000e3,y0=0e3, r_p=200e3, rho_m=3400, alpha=3e-5, g=9.8, delta_T=250, eta0=1e21, gamma=3.1
    set Function expression     = if(z==0e3, vx_bot, vx_top);0;if(z==0e3,if(((x-x0)^2+(y-y0)^2)<r_p^2, rho_m*alpha*g*r_p^2*delta_T/4/gamma/eta0*exp(-gamma*(((x-x0)^2+(y-y0)^2)/r_p^2))*365*24*60*60,0),vz_top)
  end
end

But you’re right. I should also add the lithosphere left and lithosphere right to ‘Additional tangential mesh velocity boundary’.

For the second case, I prescribed velocity on the top boundary and left the left and right boundary ‘open’ (initial lithostatic pressure) so that material can flow in from the left and out of the right boundary. The boundary condition is shown below:

subsection Boundary velocity model
set Prescribed velocity boundary indicators = bottom: function, top: function, front y: function, back y: function
  subsection Function
    set Variable names          = x,y,z
    set Function constants      = vx_top=0.02, vx_bot=0.0, vy_top=0, vy_bot=0, vz_bot=0.0, vz_top=0.0, x0=1000e3,y0=0e3, r_p=200e3, rho_m=3400, alpha=3e-5, g=9.8, delta_T=250, eta0=1e21, gamma=3.1
    set Function expression     = if(z==0e3, vx_bot, vx_top);0;if(z==0e3,if(((x-x0)^2+(y-y0)^2)<r_p^2, rho_m*alpha*g*r_p^2*delta_T/4/gamma/eta0*exp(-gamma*(((x-x0)^2+(y-y0)^2)/r_p^2))*365*24*60*60,0),vz_top)
  end
end

Yes, it is an ageing plate. The same ageing plate also applies to the first case.

subsection Initial temperature model
  set List of model names = adiabatic,function
    subsection Adiabatic
    # This is the age of the top boundary layer in years.
    set Age top boundary layer = 80e6
  end
  subsection Function
    set Variable names      = x,y,z
    set Function constants  = DeltaT=250, r_p=200e3,x0=1000e3,y0=0e3
    set Function expression = if(z>=0 && z<=340e3,if(((x-x0)^2+(y-y0)^2)<r_p^2, DeltaT*(1-((x-x0)^2+(y-y0)^2)/r_p^2),0),0)
  end
end
subsection Heating model
  set List of model names       = adiabatic heating
end

This is how the dynamic topography output looks like at 297 Myrs:


The pressure doesn’t change throughout the whole time (300 Myrs) and looks like this:

Thank you very much!

Cheers,
Ziqi