Great news! Thanks for letting us know 
A puzzling twist: I used an input thermal model with a lithosphere that goes smoothly to 273K at the surface over a distance that is well-resolved by the model mesh, set the surface temperature boundary condition to 273K, and the problem came back. But when I set the surface temperature boundary condition to 1573, even though the input thermal structure has 273 at the surface, the problem went away again. So the problem is solved when I set the bottom and top boundary temperatures to 1573, regardless of what the input thermal structure has at the boundaries.
Other than the input .txt file for the thermal structure and the Outer temperature boundary condition, the two other places I specify 1573 are the Adiabatic surface temperature to compute the reference profile, and the Reference temperature in the Visco Plastic material model subsection. I tried changing the Adiabatic surface temperature to 273, but that didn’t solve the problem, and did make the solver take a lot longer since that’s far away from the ambient mantle temperature.
Hi Sam,
Odd indeed! Can you post images of the boundary layer viscosity and temperature with the mesh also visible?
As in previous tests, does decreasing the min/max viscosity (or using a constant viscosity) remove the issue?
Cheers,
John
Hi John,
Yes of course - here are some Paraview screenshots showing the lithosphere, for the model run that did not have the “shrinking earth.” You can see the smoothly decreasing temperature towards the surface within the lithosphere, but there’s a 1-element boundary layer that is created by the 1573 K Outer temperature boundary condition. The input thermal structure has 273 K at the surface. When I set the Outer temperature to 273, the thermal structure looks as it should, with the temperature smoothly going to 273 at the surface. But this makes the “shrinking earth” come back.
Best,
Sam
Hi Sam,
Indeed, the boundary layer is certainly well resolved. In this case my guess is the issue lies within the postprocessor or how it is being used with this specific model.
That is to say unless there are other indicators (e.g., viscosity or velocity anomalies, poor solver convergence), I don’t see anything physically wrong with the initial conditions for the model with the surface temperature at 273 K.
My recollection is that many studies calculate dynamic topography using stresses from the base of the lithosphere, rather than at the surface. This would require changing where the dynamic topography post processor calculates dynamic topography, or manually calculating the values after extracting the stresses at a specific depth, isotherm, etc.
Is there a prior study that is similar in design to what you are trying to achieve here, which we could try benchmarking against?
Cheers,
John
Hi Sam,
Can you remind me whether this happens only if you use the viscoplastic rheology or also for the linear rheology? If it does happen with the linear rheology as well, we should have some runs that you could compare things to to hopefully figure our what’s going on.
Best,
Jacky
Hi Jacky,
It only happens with the Visco Plastic rheology; it works fine with the linear rheology. I’m working on making my own Python post-processor to extract the pressure at the base of the lithosphere to compute dynamic topography that way and see what results that gives.
Best,
Sam
Hi Sam,
Okay sounds good! One of my postdocs extracted the deviatoric stress tensors and the associated traction field at the base of the lithosphere using existing postprocessors to then use that as a boundary condition in a separate lithosphere deformation code (in this paper). Happy to connect you to him if that is useful at any point.
Best,
Jacky
Well it took me a while to get it all figured out, but I extracted the pressure from the base of the lithosphere, subtracted the lithostatic pressure computed from the thermal/density profile, and divided the difference by rho*g to get (air-loaded) dynamic topography. The dynamic topography field computed this way is pretty much identical to what comes out of the ASPECT post-processor (which uses surface stresses), so that’s good! I think the correspondence is because of the strong lithosphere, which transmits stresses from base to surface without much internal deformation. Unfortunately, the “shrinking earth” problem persists with this method of dynamic topography calculation as well, with the average rate of change is about the same as from the ASPECT post-processor. So this means that it’s not some quirk of the post-processor, but something more fundamental to the solution.
If I subtract the mean from the DT field at each timestep, to explicitly force a zero average rate-of-change, is that kosher? Or would I be missing something physical if I do that?
Hi Sam,
Great to hear that the new approach gives the same result at the ASPECT post processor.
Are you by chance available to join the ASPECT bi-weekly user meetings on Wednesdays (8-9 am Pacific)? The next meeting will be on November 9th. I think it would help to talk through all of this live and get input from others.
If I subtract the mean from the DT field at each timestep, to explicitly force a zero average rate-of-change, is that kosher? Or would I be missing something physical if I do that?
My initial thought is that with sufficiently stricter solver tolerances, the global rate of change should be very close to zero. Again, I think we have reached the point where a live recap with the broader group would help
Cheers,
John
Dear All,
Was there any update on this issue from the meeting? I am encountering a similar phenomenon with TERRA model output.
I was hopeful that taking stress at the base of the lithosphere would be the solution, and that the cooling lithosphere might be responsible for the net negative DT, so would like to hear if there was consensus on why that didn’t work!
Best,
Conor O’Malley
Hi Conor,
I had a baby at the beginning of November, which put a pause on my scientific progress. So I have not done any further digging since the last messages on here. But I’m intrigued that you’re having a similar issue with a different code. What sort of rheology are you using in your models?
Best,
Sam
Hi all,
First, congratulations Sam
Second, Conor thank you for posting to the forum and welcome. Despite not knowing the underlying root cause for the issue, perhaps somewhat encouraging that the phenomenon is not just specific to ASPECT.
Perhaps a path forward would be to setup a common benchmark run with ASPECT, TERRA, and CitComS and see if there is a correlation between negative DT and code/numerical method?
Cheers,
John
Hi John & Sam,
Congrats Sam!
We were observing the offset for visco-elastic runs with & without temperature-dependent viscosity. I am not running the models myself; they’re being run by Huw Davies and James Panton in Cardiff. Having spoken to them, I think the issue was we weren’t using the appropriate reference density profile to calculate normal stresses; we have resolved that issue now and get ~0 mean DT. To calculate DT, the TERRA models re-run the requested timestep instantaneously, relax plate motion constraints and impose free slip surface boundary conditions, then calculate radial normal stresses. Those are then converted to DT via simple isostatic scaling. Is that similar to how the ASPECT post-processor calculates DT?
I think benchmarking has been done before to some degree by e.g. Liu & King (2019), who benchmarked DT kernels in CitComS, ASPECT and TERRA. However, they published min/max DT and not mean values. I may have a more careful read and see if we could add anything to their results, what do you think?
Best,
Conor