[WIP] Add a cookbook that combines surface loading, viscoplastic rheologies and two-phase melting #6367
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
naliboff
requested changes
Jun 12, 2025
There was a problem hiding this comment.
@PrajaktaPMohite - Thanks for the contribution! A few initial suggestions on the PRM file.
| @@ -0,0 +1,374 @@ | |||
| # This cookbook is simple test for viscoplastic surface deformation | |||
There was a problem hiding this comment.
Suggested change
| # This cookbook is simple test for viscoplastic surface deformation | |
| # This cookbook illustrates how to combine time-dependent surface loading, free surface | |
| # deformation, a nonlinear viscoplastic rheology for the lithosphere and asthenosphere, and coupled | |
| # two-phase melting using the reactive fluid transport material model. The class of simulation | |
| # can be applied to tectonic settings where temporal changes in climate-driven surface | |
| # loading has been proposed as a mechanism to explain variations in magmatism, such | |
| # as in Clerc et al. (2024), https://doi.org/10.1038/s41467-024-45890-z. |
| @@ -0,0 +1,374 @@ | |||
| # This cookbook is simple test for viscoplastic surface deformation | |||
| # It uses a viscoplastic model to simulate surface deformation under a constant load | |||
There was a problem hiding this comment.
Suggested text for the remaining lines:
In detail, this cookbook illustrates how a time-dependent change in a parabolic surface load effects
lithospheric deformation and melting rates. This accomplished through two steps, where the first
step involves equilibration phase that lasts until the deformation and melt fields are in equilibrium
with the applied surface load. Following this phase, the surface load is decreased through time.
| set End time = 1 #1e7 | ||
| set Use years in output instead of seconds = true | ||
|
|
||
| # Melt migration models have a non linear relationship with Stokes system, temperature, |
There was a problem hiding this comment.
Suggested change
| # Melt migration models have a non linear relationship with Stokes system, temperature, | |
| # Melt migration models have a non linear relationship with Stokes system, temperature, |
| # and the advection equation for the porosity (several material properties, such as the | ||
| # viscosities and the permeability depend nonlinearly on the porosity; and changes in | ||
| # temperature determine how much material is melting or freezing). | ||
| # Here we use nonlinear solver scheme ('iterated Advection and Stokes') that iterates |
There was a problem hiding this comment.
Suggested change
| # Here we use nonlinear solver scheme ('iterated Advection and Stokes') that iterates | |
| # Here we use nonlinear solver scheme ('iterated Advection and Stokes'), which iterates |
| # viscosities and the permeability depend nonlinearly on the porosity; and changes in | ||
| # temperature determine how much material is melting or freezing). | ||
| # Here we use nonlinear solver scheme ('iterated Advection and Stokes') that iterates | ||
| # between all of these equations, and we have to set its solver tolerance and the maximum |
There was a problem hiding this comment.
Suggested change
| # between all of these equations, and we have to set its solver tolerance and the maximum | |
| # between all of these equations. We set its solver tolerance and the maximum |
| set Tangential velocity boundary indicators = left, right, bottom | ||
| end | ||
|
|
||
| # Prescribe a fixed vertical traction on the top boundary |
There was a problem hiding this comment.
Please add some documentation describing what type of load (parabolic) is being applied.
|
|
||
| ##################### Settings for melt transport ######################## | ||
|
|
||
| subsection Discretization |
There was a problem hiding this comment.
I would remove all of this, unless you are changing the default values.
|
|
||
| # The compositional fields represent: | ||
| # 1-2. fields representing porosity and peridotite | ||
| # 3-6. crustal layers |
There was a problem hiding this comment.
Suggested change
| # 3-6. crustal layers | |
| # 3-6. Layers representing rock types, which are tracked with particles. |
| end | ||
|
|
||
| # The compositional fields represent: | ||
| # 1-2. fields representing porosity and peridotite |
There was a problem hiding this comment.
Suggested change
| # 1-2. fields representing porosity and peridotite | |
| # 1-2. fields representing porosity and peridotite, which are tracked with compositional fields. |
| # Surface Heat Flow - lower crust (qs2) = 0.04533 W/m^2; | ||
| # - mantle (qs3) = 0.04033 W/m^2; | ||
| # - asthenosphere (qs4)= 0.04033 W/m^2; | ||
| # Note: The continental geotherm initial temperature model |
There was a problem hiding this comment.
Remove these lines (250-254)
naliboff
reviewed
Jun 13, 2025
There was a problem hiding this comment.
@PrajaktaPMohite - A few more suggested modifications for the PRM file
| # equations are based on the assumption that the pressure is zero at the free | ||
| # surface. | ||
| set Pressure normalization = no | ||
| set Adiabatic surface temperature = 1867.36 |
There was a problem hiding this comment.
This parameter refers to the surface adiabatic temperature, which is typically similar to the LAB temperature
| end | ||
|
|
||
| # Mesh deformation and free surface | ||
| # Advecting the free surface using a normal, rather than vertical, |
| set Initial adaptive refinement = 0 | ||
| set Initial global refinement = 5 | ||
| set Time steps between mesh refinement = 0 | ||
| set Strategy = minimum refinement function, boundary |
| # set Variable names = x,y | ||
| # set Function expression = if ( y>=150e3 && x>=140.e3 && x<=260.e3, 6, 3) | ||
| # end | ||
| # minimum of 5 global refinements |
There was a problem hiding this comment.
@PrajaktaPMohite - Don't forget to remove lines 103-111
| end | ||
|
|
||
| # Prescribe a fixed vertical traction on the top boundary | ||
| # Adding parabolic ice load on top boundary at the middle of the model. The ice load is linearly |
There was a problem hiding this comment.
Suggested change
| # Adding parabolic ice load on top boundary at the middle of the model. The ice load is linearly | |
| # Apply a parabolic load on the top boundary in the model center, which linearly |
|
|
||
| # Prescribe a fixed vertical traction on the top boundary | ||
| # Adding parabolic ice load on top boundary at the middle of the model. The ice load is linearly | ||
| # decreasing over time. |
There was a problem hiding this comment.
Suggested change
| # decreasing over time. | |
| # decreases over time. This is conceptually motivated by a melting ice sheet (e.g., Clerc et al., 2024). |
| end | ||
|
|
||
| subsection Particles | ||
| set Time between data output = 500.e3 |
There was a problem hiding this comment.
Suggested change
| set Time between data output = 500.e3 | |
| set Time between data output = 100.e3 |
| # Global parameters | ||
| set Dimension = 2 | ||
| set Start time = 0 | ||
| set End time = 1 #1e7 |
There was a problem hiding this comment.
Suggested change
| set End time = 1 #1e7 | |
| set End time = 1e6 |
PrajaktaPMohite
force-pushed
the
visco_plastic_surface_deformation_cookbook
branch
from
1a59ee5 to
e075291
Compare
PrajaktaPMohite
force-pushed
the
visco_plastic_surface_deformation_cookbook
branch
2 times, most recently
from
3fdfcd0 to
10184a0
Compare
naliboff
requested changes
Jun 14, 2025
There was a problem hiding this comment.
@PrajaktaPMohite - Thanks for adding the new text, which is in good shape. Only a few minor suggestions and comments to address. Let me know if you have any questions or have questions about the next steps.
| *This section was contributed by Prajakta Mohite & John Naliboff.* | ||
|
|
||
| This cookbook demonstrates interactions between time- and spatially-dependent surface loading,lithospheric deformation, and reactive melt transport. It is motivated by observations and prior modeling investigations highlighting the role of glacial loading and unloading on volcanic activity, and designed to serve as a potential template for investigating the effects of glacial cycles on tectonic and magmatic processes, such as in Antarctica and Iceland. | ||
| For example, the model in this cookbook can be used or extended to investigate the following topics: |
There was a problem hiding this comment.
Suggested change
| For example, the model in this cookbook can be used or extended to investigate the following topics: | |
| For example, the model in this cookbook can be used or extended to investigate the following topics: |
Add a space above this line
| This is defined in the parameter file as: | ||
| ```{literalinclude} traction_boundary.part.prm | ||
| ``` | ||
| where rho, g, and thickness represent ice density, gravitational acceleration, and ice thickness, respectively, and the spatial function defines the parabolic load geometry. Respectively, the model side and bottom boundaries contain no-slip and free-slip velocity boundary conditions. |
There was a problem hiding this comment.
Can you put rho, g, and thickness in italics (or mark them as code, can't recall exactly how to do this but there are examples in other cookbook .md files).
| ``` | ||
| where rho, g, and thickness represent ice density, gravitational acceleration, and ice thickness, respectively, and the spatial function defines the parabolic load geometry. Respectively, the model side and bottom boundaries contain no-slip and free-slip velocity boundary conditions. | ||
|
|
||
| The highly nonlinear system of equations is solved with the nonlinear solver scheme 'iterated Advection and Stokes', following prior investigations and cookbooks (see {ref}sec:cookbooks:global_melt) modeling two-phase reactive melt transport in ASPECT {cite:t}`dannberg:heister:2016`. Respectively, the linear and nonlinear solver tolerance are set to, although we note that the model dynamics may slightly or moderately differ when using stricter values. Here, the values are selected as a compromise between accuracy and model run times. We note that the introduction of a free surface, as compared to a free-slip top boundary, introduces significant challenges for both the linear and nonlinear solvers. Likewise, we note that using an open traction boundary at the model base with the traction magnitude equal to the adiabatic pressure may also introduce significant nonlinear solver issues. |
There was a problem hiding this comment.
Above, can you add in the values for the linear and nonlinear solver tolerance at the end of the sentence Respectively, the linear and nonlinear solver tolerance are set to,
|
|
||
| The highly nonlinear system of equations is solved with the nonlinear solver scheme 'iterated Advection and Stokes', following prior investigations and cookbooks (see {ref}sec:cookbooks:global_melt) modeling two-phase reactive melt transport in ASPECT {cite:t}`dannberg:heister:2016`. Respectively, the linear and nonlinear solver tolerance are set to, although we note that the model dynamics may slightly or moderately differ when using stricter values. Here, the values are selected as a compromise between accuracy and model run times. We note that the introduction of a free surface, as compared to a free-slip top boundary, introduces significant challenges for both the linear and nonlinear solvers. Likewise, we note that using an open traction boundary at the model base with the traction magnitude equal to the adiabatic pressure may also introduce significant nonlinear solver issues. | ||
|
|
||
| The initial temperature field is defined using a depth-dependent conductive temperature profile through the lithosphere in a similar fashion to the continental extension cookbook, and an adiabatic gradient of 0.5 C/km from the base of the lithosphere (80 km depth) to the model base (200 km depth). Although temperatures are specified for the model sides, these values are not used, as the boundaries are insulating (zero net heat flux). Temperature evolves in the model domain through flow induced by spatiotemporal variations in surface loading, and includes the effects of adiabatic and shear heating introduced through the use of the extended Boussinesq approximation. |
There was a problem hiding this comment.
Can you add a link to the continental extension cookbook in the sentence ... in a similar fashion to the continental extension cookbook, ...?
I think you can do this as follows:
[cookbooks/continental_extension/](https://github.com/geodynamics/aspect/tree/main/cookbooks/continental_extension)
|
|
||
| The highly nonlinear system of equations is solved with the nonlinear solver scheme 'iterated Advection and Stokes', following prior investigations and cookbooks (see {ref}sec:cookbooks:global_melt) modeling two-phase reactive melt transport in ASPECT {cite:t}`dannberg:heister:2016`. Respectively, the linear and nonlinear solver tolerance are set to, although we note that the model dynamics may slightly or moderately differ when using stricter values. Here, the values are selected as a compromise between accuracy and model run times. We note that the introduction of a free surface, as compared to a free-slip top boundary, introduces significant challenges for both the linear and nonlinear solvers. Likewise, we note that using an open traction boundary at the model base with the traction magnitude equal to the adiabatic pressure may also introduce significant nonlinear solver issues. | ||
|
|
||
| The initial temperature field is defined using a depth-dependent conductive temperature profile through the lithosphere in a similar fashion to the continental extension cookbook, and an adiabatic gradient of 0.5 C/km from the base of the lithosphere (80 km depth) to the model base (200 km depth). Although temperatures are specified for the model sides, these values are not used, as the boundaries are insulating (zero net heat flux). Temperature evolves in the model domain through flow induced by spatiotemporal variations in surface loading, and includes the effects of adiabatic and shear heating introduced through the use of the extended Boussinesq approximation. |
There was a problem hiding this comment.
Suggested change
| The initial temperature field is defined using a depth-dependent conductive temperature profile through the lithosphere in a similar fashion to the continental extension cookbook, and an adiabatic gradient of 0.5 C/km from the base of the lithosphere (80 km depth) to the model base (200 km depth). Although temperatures are specified for the model sides, these values are not used, as the boundaries are insulating (zero net heat flux). Temperature evolves in the model domain through flow induced by spatiotemporal variations in surface loading, and includes the effects of adiabatic and shear heating introduced through the use of the extended Boussinesq approximation. | |
| The initial temperature field is defined using a depth-dependent conductive temperature profile through the lithosphere in a similar fashion to the continental extension cookbook, and an adiabatic gradient of 0.5 C/km from the base of the lithosphere (80 km depth) to the model base (200 km depth). Although temperatures are specified for the model sides, these values are not used, as the boundaries are insulating (zero net heat flux). Temperature evolves in the model domain through flow induced by spatiotemporal variations in surface loading, and includes the effects of adiabatic and shear heating introduced through the use of the extended Boussinesq approximation. |
PrajaktaPMohite
force-pushed
the
visco_plastic_surface_deformation_cookbook
branch
4 times, most recently
from
438a2a8 to
f9237db
Compare
PrajaktaPMohite
force-pushed
the
visco_plastic_surface_deformation_cookbook
branch
from
f9237db to
74c4e9e
Compare
xlia62
reviewed
Jun 20, 2025
| # diffusion or composite flow law. Here, dislocation is selected | ||
| # so no need to specify diffusion creep parameters below, which are | ||
| # only used if "diffusion" or "composite" option is selected. | ||
| set Viscous flow law = dislocation |
There was a problem hiding this comment.
viscous flow law has already been defined beofre as composite
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This cookbook provides an example for how to combine surface loading, nonlinear lithospheric and asthenospheric deformation, and two-phase reactive melt transport using the reactive fluid transport material model.