Hi all,
I’m using quasi-instantaneous ASPECT models (a few short timesteps) to estimate the rate of change of dynamic topography in a spherical Earth-like model. However, I’m running into an unexpected problem. When I use the visco-plastic rheology (which is what I would like to use for a realistic model), I find that dDT/dt is negative over the large majority of the surface, and the spherically averaged value is negative. I don’t find this when using the simple rheology, which gives an average dDT/dt of 0, or with the simple temperature-dependent rheology. I get the non-zero average with diffusion creep only as well as with diffusion+dislocation, so it’s not just the non-linearity of dislocation creep that’s causing this.
I’m not sure if this is a numerical issue or a physics problem. Obviously the Earth isn’t shrinking (other than very slowly due to heat loss), but could this be the result of computing DT from stresses? Should we expect that the surface normal stress should average to 0, or that the average shouldn’t change over time? My intuition is yes, at least for a linear rheology, but I’m not positive. Alternatively, if it’s a numerical artifact, is there a simple fix?
Best,
Sam
1 Like
Hi Sam,
Thanks for posting this to the forum. It sounds like the issue is related to the addition of plasticity.
Dynamic topography is calculated using the dynamic pressure (i.e., full pressure - adiabatic pressure profile). The values are calculated everywhere in ASPECT, but what the post processor extracts is the values at the surface calculated from the topmost element.
If negative pressures are present in the top few elements, this may lead to unrealistic dynamic pressure gradients being supplied to the dynamic topography post processor.
Here are a few ways to check if this is indeed the issue and potential fixes:
- Check the pressure field in each of the model runs to see if there are any negative values (should only be near the surface). The negative pressures may be present in all models, but only producing the observed when plasticity is introduced.
- Set
Pressure normalization = surface and then set Surface pressure equal to a nonzero value. This will simply scale all the pressures in the model.
- Make sure
Allow negative pressures in plasticity = false. This will supply the Drucker Prager yield criterion a pressure of 0 if it is negative.
Cheers,
John
Hi John,
Thanks for your quick response! I don’t want plasticity in my model, so I have tried to disable it by setting the cohesion to a very high value (1e20). But it could still be throwing it off. I will try your proposed fixes and see how it goes.
Best,
Sam
Hi Sam,
Ah, apologies, I misread your first post!
Indeed, if the cohesion is set to 1e20, then it is likely dislocation creep producing the odd values. However, I still suspect negative pressures might be the culprit here. If that is the case, the fix on point 2 (pressure normalization) may fix the issue.
Cheers,
John
Hi John,
I ran a model run where I set the average surface pressure to 1.5e9, which is well in excess of the largest negative pressures I was getting with it set to 0, and successfully eliminated negative pressures from the near-surface. Unfortunately, it didn’t fix the shrinking problem - there’s still an average decrease in dynamic topography of ~0.2 mm/yr, similar to the same run with a surface pressure set to 0, and 95% of the surface has a negative dDT/dt.
Best,
Sam
Hi Sam,
Thanks for running that test, and interesting that the issue is still persisting.
Question 1: Just to be sure, the model is incompressible and you are using the Visco Plastic material model?
Question 2: Re-reading your first post, it sounds like you have a model with diffusion + dislocation creep that gives an average dynamic topography of 0 through time, and the current behavior only arises when you add plasticity (even though the cohesion is set to a high value). Is this correct? If not, can specify exactly what changes to the .prm file produce the issue?
Question 3 - Is the model converging to reasonable linear and non-linear solver tolerance values? Would you mind posting the details of the solver parameters? I don’t think this is a likely culprit, but it can’t hurt to double check?
Cheers,
John
Hi John,
Thanks for continuing to help me troubleshoot. I am using the incompressible “Visco Plastic” material model when I get the shrinking earth. I have also tried the “Ascii reference profile” material model (with and without lateral viscosity variations) and the “Simple” material model, without getting the problem.
I have used the Visco Plastic model for diffusion only (by setting the prefactor for dislocation to be very low and the cohesion to be very high) as well as diffusion+dislocation (by setting the cohesion to be very high); both scenarios cause the problem. I have not tried enabling plasticity by using reasonable parameters for that part. To summarize, I have gotten the shrinking earth with every run where I’ve used “Visco Plastic”, but not with simpler rheologies. Below you can see the material model parameters I’m using for the diffusion+dislocation scenario.
As regards the solver parameters, I’ve pasted them below as well. I don’t really know if they’re reasonable values or not, but I would love a more expert opinion. I did have to relax them as I added complexity to the model because it wasn’t converging, but I don’t know if that’s a result of the material model or of the short-wavelength mantle structure I’m trying to model.
Thanks again,
Sam
### SOLVER PARAMETERS ###
# non-linear solver
set Nonlinear solver scheme = single Advection, iterated Stokes
set Nonlinear solver tolerance = 1.0e-3
set Max nonlinear iterations = 500
# linear solver
subsection Solver parameters
subsection Stokes solver parameters
set Linear solver A block tolerance = 1e-1
set Linear solver tolerance = 1e-3
set Maximum number of expensive Stokes solver steps = 4000
set Number of cheap Stokes solver steps = 300
set GMRES solver restart length = 200
end
set Temperature solver tolerance = 1e-7
end
subsection Material model
set Model name = visco plastic
subsection Visco Plastic
# reference stuff
set Reference temperature = 1573
set Minimum strain rate = 1.e-20
set Adiabat temperature gradient for viscosity = 9.24e-09 # 0.3 K/km
# Set up phase transitions so that dislocation becomes inactive in the lower mantle for all fields.
set Phase transition depths = background: 660e3
set Phase transition widths = background: 1e3
set Phase transition temperatures = background: 1573
set Phase transition Clapeyron slopes = background: 0
# Viscosity cut offs
set Minimum viscosity = background: 1.0e19|1.0e20
set Maximum viscosity = background: 1.0e23|1.0e24
# density-related parameters
set Thermal diffusivities = 1.e-6
set Heat capacities = 750
set Densities = 3300
set Thermal expansivities = 3e-5
# viscosity-related parameters
set Viscous flow law = composite
set Viscosity averaging scheme = harmonic
# dislocation creep (just z < 660 km)
# We set the prefactor to a very low value in the lower mantle to effectively turn off dislocation creep.
set Prefactors for dislocation creep = background: 1.640436e-15|1e-40
set Stress exponents for dislocation creep = 3.5
set Activation energies for dislocation creep = 540.e3
set Activation volumes for dislocation creep = 12e-6
# diffusion creep, (gives lower/upper mantle viscosity = 10)
set Prefactors for diffusion creep = background: 6.375796e-11|7.231448e-13
set Stress exponents for diffusion creep = 1
set Grain size exponents for diffusion creep = 0
set Activation energies for diffusion creep = 300e3
set Activation volumes for diffusion creep = background: 4e-6|2.5e-6
# Plasticity parameters using drucker-prager - turn off using high cohesion
set Cohesions = background: 1.0e20|1.0e20
end
end
Hi Sam,
Thanks for providing all of that information and hopefully we are bit closer to debugging the issue.
Your choice of rheology parameters looks reasonable, but I do have some suggestions for modifying the solver settings.
Here are some ideas for what to try next:
-
Try setting the min/max viscosity to the same value (ex: 1.e22 Pa s), so that you have a uniform viscosity everywhere. If that produces a reasonable spatial and time-dependent patterns of dynamic topography, try systematically increasing the min/max viscosity contrast by one order of magnitude to see how this affects the patterns.
-
After the tests in point 1, try using the solver parameters listed below. At some point it would be good to switch to using block GMG, but I need to check with the other developers on what tweaks need to be made for spherical models (null space removal, etc).
Cheers,
John
set Nonlinear solver scheme = single Advection, iterated defect correction Stokes
subsection Solver parameters
subsection Stokes solver parameters
set Stokes solver type = block AMG
set Number of cheap Stokes solver steps = 1000
set Use full A block as preconditioner = true
set Linear solver tolerance = 1e-7
end
subsection Newton solver parameters
set SPD safety factor = 0.9
set Maximum linear Stokes solver tolerance = 1e-2
set Use Eisenstat Walker method for Picard iterations = true
end
end
Hi John,
Thanks for the ideas. I submitted a set of 5 runs where the allowable range of viscosity spans 0, 1, 2, 3, and 4 orders of magnitude. I’ll update you with the results of that when I have them, then try changing the solver parameters.
Best,
Sam
Hi John,
For the run with no range of viscosity (actually a min of 1e21 and a max of 1.0001e21), I got a very slight negative average dDT/dt (~1e-5 m/yr). For 1 order of magnitude it was ~4e-5, and for two it was ~8e-5. It didn’t really increase beyond that (this is using a simpler mantle structure than the 2e-4 m/yr I referenced above). I’ll try your suggested solver parameters now and see what happens.
Best,
Sam
Hi Sam,
Thanks for the update. Encouraging to know that there is some correlation between the dDT/dt patterns and rheology.
A few other ideas to consider, depending on the outcome of the next tests:
- Just to be sure, is the density formulation used in the
Visco Plastic model producing the same density field as with other material models?
- Can you modify the flow law equations parameters for the
Visco Plastic model to produce the same rheology structure used in the other material models? For example, removing strain-rate dependence or even using a constant viscosity (set n=1, E=0, V=0, etc)?
Cheers,
John
Hi John,
Sorry for the delayed response. Things got hung up in the cluster queue. Unfortunately, I’m not getting the solver to converge with the parameters you suggested. Error message:
The iterative Stokes solver in Simulator::solve_stokes did not
converge.
The initial residual was: 1.722873e+20
The final residual is: 5.256358e+15
The required residual for convergence is: 1.722873e+13
Iterative method reported convergence failure in step 1000. The
residual in the last step was 5.25636e+15.
I tried again allowing 5000 steps, and got the same thing, with essentially the same residual.
The iterative Stokes solver in Simulator::solve_stokes did not
converge.
The initial residual was: 1.722873e+20
The final residual is: 5.174458e+15
The required residual for convergence is: 1.722873e+13
Iterative method reported convergence failure in step 5000. The
residual in the last step was 5.17446e+15.
Here’s the solver chunk of the parameter file that I was using:
set Nonlinear solver scheme = single Advection, iterated defect correction Stokes
subsection Solver parameters
subsection Stokes solver parameters
set Stokes solver type = block AMG
set Number of cheap Stokes solver steps = 2000
set Use full A block as preconditioner = true
set Linear solver tolerance = 1e-7
set Maximum number of expensive Stokes solver steps = 5000
end
subsection Newton solver parameters
set SPD safety factor = 0.9
set Maximum linear Stokes solver tolerance = 1e-2
set Use Eisenstat Walker method for Picard iterations = true
end
end
I’m using a simple tomography scaling for the thermal (and thus viscosity) structure, so there shouldn’t be any weird singularities or unusually high gradients in parameters. Should I relax the tolerances, or would that result in an unacceptably inaccurate solution?
Best,
Sam
Hi Sam,
Can you post the entire contents of the log.txt file? I’m curious to see the convergence behavior throughout the model run.
Going forward, I think it may be easiest to take a step back and try to isolate the issue in 2D models (cartesian or spherical).
From your current tests, it seems the issue is with the material model and thus hopefully the issues would be present in 2D simulations as well.
This will allow for faster debugging and allow me (and others) to easily run the tests as well.
If this sounds reasonable, can you come up with a 2D model to test the different material models on? I think starting with cartesian would be best to further simplify the problem. Perhaps even something like the one of the simple convection box models would work? They all should produce dT/dt = 0 if the equations are being solved correctly.
Thanks for your continued patience through the debugging process 
Cheers,
John
Hi John,
Here’s the log and the solver history for the last run (with 5000 steps).
solver_history.txt (114.0 KB)
log.txt (16.8 KB)
I’ll work on replicating the problem in a 2-D model and update you accordingly.
Best,
Sam
Hi John,
I think I’ve more or less replicated the problem in 2-D. I tried to keep it similar to the 3-D spherical model I was doing, but I simplified the mantle structure to two antipodal anomalies, one negative, one positive (see here for input text file). I’ve done two rheologies, one “simple” (uniform), and one with “Visco Plastic” (diffusion and dislocation creep). For each I did both your suggested solver parameters and the ones I was using before, although there was no major difference in the results. Parameter files for each are pasted below.
Both rheologies produce plausible flow field, namely upwelling for the hot/buoyant anomaly and downwelling for the cold/dense anomaly. The dynamic topography field corresponds, with higher DT over the upwelling and lower DT over the downwelling. Curiously, the DT field does not average to 0 in either case. With a simple rheology, the average is -440 m, with peaks at -1790/+912. With the visco plastic rheology, the average is -816 m, with peaks at -1130/-437 (i.e. DT is negative everywhere). Other than the offset mean, the DT fields look symmetric, as expected.
As regards the rate of change: For the simple rheology, I find an average dDT/dt of ~1e-6 m/yr (net uplift). For the visco-plastic rheology, I get an average dDT/dt of ~-2.5e-5 m/yr, so an order of magnitude more, and change in sign to be net subsidence/shrinking. The good news is that these runs take 1-2 minutes on a single node, so it’s easy for me (or anyone else) to try things quickly and easily.
Best,
Sam
annulus_simple_Glob3Adapt2.prm (4.7 KB)
annulus_simple_Glob3Adapt2_newsolver.prm (4.7 KB)
annulus_disl_Glob3Adapt2.prm (6.2 KB)
annulus_disl_Glob3Adapt2_newsolver.prm (6.3 KB)
Here’s what the fields look like for dynamic topography and its rate of change:
Hi Sam, hi John,
I think the problem, at least with the 2D model, is with the boundary temperature model. The outer and inner boundaries are fixed to be 273 and 2273 K, but both boundaries have a temperature of 1573 K according to the initial temperature model.
This isn’t the first time we’ve seen poor behaviour when the boundary conditions are inconsistent with the conditions inside the domain. It might be instructive to see how quickly the problem goes away if you smoothed the temperature jump over a few cells.
Best wishes,
Bob
Hi Bob,
Thanks for your insights. I just did a run where I set both of the boundary conditions to 1573 (the same as the ambient mantle), and lo and behold the problem went away! It gives an average DT of -0.07 m and average dDT/dt of 1e-8 m/yr. Since I’m interested in nearly-instantaneous models, I don’t really need to worry about heat flow across the boundaries. I’ll try slowly adding back levels of complexity and see if this fix still completely solves the problem. Thanks!
Best,
Sam
Hi Sam, Hi Bob,
Great news, thanks for the sharp eye Bob!
Sam - please let us know how the additional tests go with adding more complexity back in, especially if anything seems off in the behavior as the rheology gets more complicated.
Cheers,
John
Hi John and Bob,
Good news! I did a 3-D run using a simple tomography scaling for the mantle structure, with revised boundary conditions to be 1573 at both top and bottom, which is the reference temperature I used to scale the input thermal structure, and it seems to be behaving well. This is with the visco-plastic rheology using diffusion and dislocation creep.
My next step is to try my more complex thermal input model. This one differs in that I impose a lithosphere with a half-space cooling profile that’s 273K at the surface, so the 273 surface boundary condition is still appropriate. I don’t have any increase in temperature at the bottom boundary condition, so I will change the 2273 to 1573 and cross my fingers.
Best,
Sam
