New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add option to change viscosity and/or density of composition spcrust #2586
Conversation
/run-tests |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Partial review through the visco_plastic.h and visco_plastic.cc. Mostly minor comments regarding documentation and variable names.
As we discussed offline, this is a fairly specific feature that adds a fair bit of code. While I think it will be to a fair number of people, I am still a bit concerned that if it is included it should be expanded to work on more than one compositional field. For example, what about models where there are multiple layers representing crustal material.
I'd would like to get the other developers (and users) thoughts on this before proceeding further.
/** | ||
* Transition from maximum spcrust viscosity to flow-law defined value | ||
* over specified pressure range | ||
*/ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
-
Fix formatting (spacing of comments)
-
Can you expand the comment a bit to make it clear that
spcrust
is the name of a compositional field that will have it's values altered over a pressure range? It may also be worth specifying exactly what each of thespcrust
variables represent below, as the namespcrust
is not particularly descriptive. -
What is the origin of the name
spcrust
? Specifically, what do the "s" and "p" stand for? For clarity, it might be worth expanding the name to include what "s" and "p" represent.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
spcrust stands for subducting-plate crust. I've expanded the comments to clarify the parameters and what the options do.
@@ -248,6 +248,28 @@ namespace aspect | |||
// Apply strain weakening of the viscous viscosity | |||
viscosity_pre_yield *= viscous_weakening; | |||
|
|||
// For composition spcrust, change from fixed *maximum* viscosity to flow-law viscosity |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Potential rewording: For composition spcrust
-> If there is a compositional field named spcrust, change the viscosity from the fixed ....
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've reworded this starting from your suggestion
@@ -391,10 +413,26 @@ namespace aspect | |||
double density = 0.0; | |||
for (unsigned int j=0; j < volume_fractions.size(); ++j) | |||
{ | |||
double delta_crust_density = 0.0; | |||
if (use_spcrust_density_change) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can you add a comment similar to the one for the spcrust
viscosity above?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
done.
else if ((pressure > spcrust_viscosity_minimum_pressure) | ||
&& (pressure < spcrust_viscosity_maximum_pressure)) | ||
{ | ||
viscosity_pre_yield = std::min(maximum_spcrust_viscosity* |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think it would be helpful to have a comment briefly writing out and explaining the formula below (it was not intuitive at first glance).
Also, comparing to the formula for the density change, is there no need for a condition when pressure >= spcrust_viscosity_maximum_pressure
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've written a full comment here to explain the calculation. We want to cap the maximum value to the spcrust viscosity at shallow depth and it linearly increases with depth, but we want to allow this crustal value to be lower if the composite rheology predicts that it is lower. There is no need for a condition from pressure above the maximum pressure because in this case the flow-law-defined value (viscosity_pre_yield) is used as is.
"from a constant value to the value determined by the flow law parameters. " | ||
"Units: None"); | ||
prm.declare_entry ("Maximum spcrust viscosity", "1e28", Patterns::Double(0), | ||
"Maximum viscosity for the composition called spcrust. Using a value of 1e20 $Pa \\, s$" |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Missing space after $
"the flow law. A value of 2.0e9 $Pa$ would correspond to a " | ||
"depth of about 60 km. Units: $Pa$"); | ||
prm.declare_entry ("Maximum transition pressure spcrust viscosity", "0.0", Patterns::Double(0), | ||
"Pressure at which to end smooth transition from " |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
"to end the smooth transition .."
"the flow law. A value of 3.9e9 $Pa$ would correspond to a " | ||
"depth of about 120 km. Units: $Pa$"); | ||
|
||
// Transition the spcrust density from defined value by a delta-rho |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
What is the "defined value" for density? Is this just the regular density defined by the reference density, temperature and thermal expansivity?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
defined value is that given by the equation of state in the material model. I've changed this to say equation of stated defined value
"For basalt density of 3000 and eclogite density of 3540, use a value of 540 $kg/m^3" | ||
"Units: $kg/m^3$"); | ||
prm.declare_entry ("Minimum transition pressure spcrust density", "0.0", Patterns::Double(0), | ||
"Pressure at which to start the smooth transition from in density." |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Missing or extra words here?
I agree with John that this might be outside of the scope of this general material model. Magali, would you need these features if you were to make a general subduction cookbook, or is it just something you or the CIDER group would need for a specific model setup? If it is the second I can show you how to use the visco plastic model as a base model for a model-specific plugin that you can keep outside of the main repository. If it is the first we should talk about if it should go into the material model, or if we should make a specific material model in the cookbooks folder. |
I've gone through and expanded the comments as requested by @naliboff. As to the question of whether this belongs in the material model or not... I would argue that these two options are required to do the most basic subduction model. It is possible to make a subduction model in which the crustal viscosity remains weak, but the dynamics of the slab are completely different than allowing this value to increase. Being able to modify how these parameters change from the parameter file allows a whole range of hypothesis and behaviors to be tested. Similarly, for the density, one could model subduction with no crustal density in the subducting plate crust, but if you include a crustal density, then you really need to include the change to eclogite since this is almost 500 kg/m^3!! Perhaps @naliboff and @gassmoeller can discuss this further and determine if this is okay as part of this material model (perhaps with the understanding that future planned changes to the material model will allow combinations of rheology/density behaviors to be put together as separate modules more easily). |
@mibillen - Thanks for adding the additional documentation! @gassmoeller and I (and possibly others) are going to discuss the proposed functionality and various options early next week. I think it might be useful to expand the proposed functionality to allow the same modifications for any compositional field representing a rock type, but I'll think about the code structure a bit more before expanding on this. More soon and thanks for the patience in waiting for us to talk this over! |
Are you thinking that instead of connecting this to just spcrust, the user could give a list of composition names, and then a matching lists for drho, minimum_pressure and maximum_pressure and then similarly matching lists of maximum_viscosity, minimum pressure and maximum pressure? |
@mibillen - Yes, this is exactly what I had in mind. With this general implementation, the change can be applied to a multi-layer ocean crust, sediment layers, etc. There may be also be case where some of these changes would be helpful even without subduction (ex: lithospheric delamination). Thoughts everyone else? |
@naliboff and @gassmoeller - I think I should close this pull request because I will use the new pressure-transition functionality to make these changes. Please confirm that this is what I should do and I will close this pull request |
Closing this pull request because we will make these changes using the phase transition capabilities being added with #2656. |
Add options to visco_plastic module for allowing the viscosity and/or density of a specific composition (must be called spcrust) to change linearly from one value to another over a specific pressure range. This can be used to model things like the change from basalt to eclogite in the subducting plate crust during subduction, but it could also be used for modeling other layers (like a serpentinite layer). Since this change doesn't conflict with the other pull request for visco_plastic, I decided to go ahead and get this one submitted also. (After this, just one more change to allow the viscosity in the lower mantle to be different than the upper mantle).