Skip to content

Track free water from each lithology in the Tian 2019 Reaction Model #6556

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.

Sign up for GitHub

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

Open
wants to merge 2 commits into
base: main
Choose a base branch
from
Open

Track free water from each lithology in the Tian 2019 Reaction Model #6556

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c8caa8a
Track free water from each lithology
danieldouglas92 Jun 18, 2025
b9e8686
add test case
danieldouglas92 Jun 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions include/aspect/material_model/reaction_model/tian2019_solubility.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,12 @@ namespace aspect
void
parse_parameters (ParameterHandler &prm);

/**
* A boolean that indicates whether the model tracks the contribution of each
* lithology to the total volume fraction of free water.
*/
bool track_unique_water_components;

private:

/**
Expand Down
25 changes: 25 additions & 0 deletions source/material_model/reaction_model/tian2019_solubility.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,6 +147,13 @@ namespace aspect
void
Tian2019Solubility<dim>::declare_parameters (ParameterHandler &prm)
{
prm.declare_entry ("Track source of free water", "false",
Patterns::Bool(),
"Whether to track the volume fraction of the free water that comes from each of the "
"four rock compositions (sediment, MORB, gabbro, peridotite). If false, only the "
"total water is tracked. If true, the total water, and the volume fraction from "
"each rock composition is tracked. This requires defining four compositional fields: "
"sediment_porosity, MORB_porosity, gabbro_porosity, and peridotite_porosity. ");
prm.declare_entry ("Maximum weight percent water in sediment", "3",
Patterns::Double (0),
"The maximum allowed weight percent that the sediment composition can hold.");
Expand Down Expand Up @@ -178,6 +185,24 @@ namespace aspect
AssertThrow(this->introspection().compositional_name_exists("peridotite"),
ExcMessage("The Tian approximation only works "
"if there is a compositional field called peridotite."));

track_unique_water_components = prm.get_bool ("Track source of free water");
if (track_unique_water_components)
{
AssertThrow(this->introspection().compositional_name_exists("sediment_porosity"),
ExcMessage("Tracking the contribution of free water sourced from sediment "
"requires a compositional field named sediment_porosity."));
AssertThrow(this->introspection().compositional_name_exists("MORB_porosity"),
ExcMessage("Tracking the contribution of free water sourced from MORB "
"requires a compositional field named MORB_porosity."));
AssertThrow(this->introspection().compositional_name_exists("gabbro_porosity"),
ExcMessage("Tracking the contribution of free water sourced from gabbro "
"requires a compositional field named gabbro_porosity."));
AssertThrow(this->introspection().compositional_name_exists("peridotite_porosity"),
ExcMessage("Tracking the contribution of free water sourced from peridotite "
"requires a compositional field named peridotite_porosity."));
}

tian_max_peridotite_water = prm.get_double ("Maximum weight percent water in peridotite");
tian_max_gabbro_water = prm.get_double ("Maximum weight percent water in gabbro");
tian_max_MORB_water = prm.get_double ("Maximum weight percent water in MORB");
Expand Down
55 changes: 52 additions & 3 deletions source/material_model/reactive_fluid_transport.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -188,7 +188,7 @@ namespace aspect

if (fluid_solid_reaction_scheme != katz2003)
{
const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");
const unsigned int porosity_idx = this->introspection().compositional_index_for_name("porosity");

// Modify the viscosity from the base model based on the presence of fluid.
if (in.requests_property(MaterialProperties::viscosity))
Expand Down Expand Up @@ -252,6 +252,32 @@ namespace aspect
std::vector<double> eq_free_fluid_fractions(out.n_evaluation_points());
melt_fractions(in, eq_free_fluid_fractions, &out);

// Initialize the indices of compositional fields as invalid ints. These will
// only be used if the Tian 2019 reaction model with tracking unique water components
// is defined in the parameter file.
unsigned int sediment_porosity_idx = numbers::invalid_unsigned_int;
unsigned int MORB_porosity_idx = numbers::invalid_unsigned_int;
unsigned int gabbro_porosity_idx = numbers::invalid_unsigned_int;
unsigned int peridotite_porosity_idx = numbers::invalid_unsigned_int;

unsigned int sediment_idx = numbers::invalid_unsigned_int;
unsigned int MORB_idx = numbers::invalid_unsigned_int;
unsigned int gabbro_idx = numbers::invalid_unsigned_int;
unsigned int peridotite_idx = numbers::invalid_unsigned_int;

if (tian2019_model.track_unique_water_components)
{
sediment_porosity_idx = this->introspection().compositional_index_for_name("sediment_porosity");
MORB_porosity_idx = this->introspection().compositional_index_for_name("MORB_porosity");
gabbro_porosity_idx = this->introspection().compositional_index_for_name("gabbro_porosity");
peridotite_porosity_idx = this->introspection().compositional_index_for_name("peridotite_porosity");

sediment_idx = this->introspection().compositional_index_for_name("sediment");
MORB_idx = this->introspection().compositional_index_for_name("MORB");
gabbro_idx = this->introspection().compositional_index_for_name("gabbro");
peridotite_idx = this->introspection().compositional_index_for_name("peridotite");
}

for (unsigned int q=0; q<out.n_evaluation_points(); ++q)
{
const unsigned int bound_fluid_idx = this->introspection().compositional_index_for_name("bound_fluid");
Expand All @@ -274,7 +300,7 @@ namespace aspect
// mass fraction of water in the solid, Fl_wt is the wt% of water in the fluid, and Fl_mass is
// the mass fraction of water in the fluid. Fl_wt is always 100%, because the fluid is assumed to
// be composed of only water. The bound_fluid composition gives Sm_wt, and we compute Sm_mass as
// was done when defining mass_frac_porosity.
// was done when defining mass_fraction_porosity.
const double total_mass_fraction_water = in.composition[q][bound_fluid_idx] * (1 - volume_fraction_porosity) * solid_density / bulk_density +
mass_fraction_porosity;

Expand All @@ -286,7 +312,7 @@ namespace aspect
// Since porosity is a volume fraction, convert the mass fraction change to a volume fraction change
// to update the porosity value. We cannot use the bulk density to do this, because we do not yet know
// the new bulk density of the rock after the reactions have been applied. We do this by substituting
// the new bulk density as a function of the new volume fraction porosity into this equation:
// the new bulk density as a function of the new volume fraction porosity into this equation:
// new_Fl_mass * new_bulk_density = new_volume_frac * fluid_density
// and solving for the new volume fraction. This gives us the following equation:
const double new_volume_fraction_porosity = solid_density * new_mass_fraction_porosity /
Expand Down Expand Up @@ -324,6 +350,29 @@ namespace aspect
reaction_rate_out->reaction_rates[q][c] = porosity_change / reaction_time_step_size;
else
reaction_rate_out->reaction_rates[q][c] = 0.0;

if (tian2019_model.track_unique_water_components)
{
// Determine what the fraction of each lithology is at the current quadrature point.
const double sediment_percent = std::min(std::max(in.composition[q][sediment_idx], 0.0), 1.0);
const double MORB_percent = std::min(std::max(in.composition[q][MORB_idx], 0.0), 1.0);
const double gabbro_percent = std::min(std::max(in.composition[q][gabbro_idx], 0.0), 1.0);
const double peridotite_percent = std::min(std::max(in.composition[q][peridotite_idx], 0.0), 1.0);

// The total porosity change is determined based on the average lithology at each point.
// Therefore to track the porosity change for each lithology, we need to multiply the
// total by the contribution of each lithology.
const double porosity_change_for_tracking = porosity_change < 0 ? 0 : porosity_change / reaction_time_step_size;

if (c == sediment_porosity_idx)
reaction_rate_out->reaction_rates[q][sediment_porosity_idx] = porosity_change_for_tracking * sediment_percent;
else if (c == MORB_porosity_idx)
reaction_rate_out->reaction_rates[q][MORB_porosity_idx] = porosity_change_for_tracking * MORB_percent;
else if (c == gabbro_porosity_idx)
reaction_rate_out->reaction_rates[q][gabbro_porosity_idx] = porosity_change_for_tracking * gabbro_percent;
else if (c == peridotite_porosity_idx)
reaction_rate_out->reaction_rates[q][peridotite_porosity_idx] = porosity_change_for_tracking * peridotite_percent;
}
}
}
}
Expand Down
79 changes: 79 additions & 0 deletions tests/tian_track_fluid_sources.prm
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,79 @@
# This model generates a representation of a phase diagram for the equilibrium
# wt% water that can be stored in an upper mantle peridotite (PUM) as shown in
# Figure 11 of Tian et. al., 2019. This model reproduces that phase diagram
# (very low res by default), with an important caveat being that the pressure axis
# is inverted relative to the ASPECT output. A surface pressure of 600 MPa is imposed
# to reproduce the pressure axis,and a lateral temperature gradient is imposed to
# reproduce the temperature axis.
include $ASPECT_SOURCE_DIR/tests/tian_peridotite.prm
set Nonlinear solver scheme = iterated Advection and Stokes
# set Surface pressure = 0.6e9

# The Tian approximation implementation requires that all 4 rock compositions exist
# as compositional fields, though they do not have to all be non zero.
subsection Compositional fields
set Number of fields = 10
set Names of fields = porosity, bound_fluid, peridotite, gabbro, MORB, sediment, peridotite_porosity, gabbro_porosity, MORB_porosity, sediment_porosity
set Compositional field methods = darcy field, field, field, field, field, field, darcy field, darcy field, darcy field, darcy field
set Types of fields = porosity, generic, chemical composition, chemical composition, chemical composition, chemical composition, porosity, porosity, porosity, porosity
end

# Specify that the entire domain is peridotite, with an initial bound water content
# of 11% everywhere, the bound_water will evolve to the parametrized phase diagram
# for peridotite.
subsection Initial composition model
set Model name = function

subsection Function
set Function constants = initial_porosity=0.0, initial_bound_water=0.11
set Function expression = initial_porosity; \
initial_bound_water; \
0.4; \
0.3; \
0.2; \
0.1; \
0; \
0; \
0; \
0
end
end

# No gravity
subsection Gravity model
set Model name = vertical

subsection Vertical
set Magnitude = 0.0
end
end

# Material model, using the simpler base model and selecting the tian approximation
# as the fluid solid reaction scheme for the reactive fluid transport model.
subsection Material model
set Model name = reactive fluid transport

subsection Reactive Fluid Transport Model
set Base model = visco plastic
set Reference fluid density = 2000
subsection Tian 2019 model
set Track source of free water = true
set Maximum weight percent water in peridotite = 11
set Maximum weight percent water in gabbro = 11
set Maximum weight percent water in MORB = 11
set Maximum weight percent water in sediment = 11
end
end

subsection Visco Plastic
set Viscous flow law = diffusion
set Densities = 3000, 3000, 3000, 3300, 3100, 2900, 2300, 3000, 3000, 3000, 3000
set Thermal expansivities = 0
set Minimum viscosity = 1e21
set Maximum viscosity = 1e21
end
end

subsection Melt settings
set Include melt transport = false
end
93 changes: 93 additions & 0 deletions tests/tian_track_fluid_sources/screen-output
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
-----------------------------------------------------------------------------
-- This is ASPECT --
-- The Advanced Solver for Planetary Evolution, Convection, and Tectonics. --
-----------------------------------------------------------------------------
-- . version 3.1.0-pre (track_fluid_source, 986c9a918)
-- . using deal.II 9.6.2
-- . with 32 bit indices
-- . with vectorization level 1 (SSE2, 2 doubles, 128 bits)
-- . using Trilinos 14.4.0
-- . using p4est 2.3.6
-- . using Geodynamic World Builder 1.0.0
-- . running in DEBUG mode
-- . running with 1 MPI process
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
-- For information on how to cite ASPECT, see:
-- https://aspect.geodynamics.org/citing.html?ver=3.1.0-pre&mf=1&sha=986c9a918&src=code
-----------------------------------------------------------------------------
Number of active cells: 16 (on 3 levels)
Number of degrees of freedom: 1,078 (162+25+81+81+81+81+81+81+81+81+81+81+81)

*** Timestep 0: t=0 years, dt=0 years
Solving temperature system... 0 iterations.
Skipping porosity composition solve because RHS is zero.
Solving bound_fluid system ... 0 iterations.
Solving peridotite system ... 0 iterations.
Solving gabbro system ... 0 iterations.
Solving MORB system ... 0 iterations.
Solving sediment system ... 0 iterations.
Skipping peridotite_porosity composition solve because RHS is zero.
Skipping gabbro_porosity composition solve because RHS is zero.
Skipping MORB_porosity composition solve because RHS is zero.
Skipping sediment_porosity composition solve because RHS is zero.
Solving Stokes system (GMG)... 0+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.41953e-16, 0, 1.14893e-16, 1.43204e-16, 8.69635e-17, 1.43204e-16, 1.43204e-16, 0, 0, 0, 0, 0
Relative nonlinear residual (total system) after nonlinear iteration 1: 1.43204e-16


Postprocessing:
Compositions min/max/mass: 0/0/0 // 0.11/0.11/4.4e+09 // 0.4/0.4/1.6e+10 // 0.3/0.3/1.2e+10 // 0.2/0.2/8e+09 // 0.1/0.1/4e+09 // 0/0/0 // 0/0/0 // 0/0/0 // 0/0/0

*** Timestep 1: t=1000 years, dt=1000 years
Solving composition reactions... in 100 substep(s).
Solving temperature system... 6 iterations.
Solving porosity system ... 0 iterations.
Solving bound_fluid system ... 0 iterations.
Solving peridotite system ... 0 iterations.
Solving gabbro system ... 0 iterations.
Solving MORB system ... 0 iterations.
Solving sediment system ... 0 iterations.
Solving peridotite_porosity system ... 0 iterations.
Solving gabbro_porosity system ... 0 iterations.
Solving MORB_porosity system ... 0 iterations.
Solving sediment_porosity system ... 0 iterations.
Solving Stokes system (GMG)... 0+0 iterations.
Relative nonlinear residuals (temperature, compositional fields, Stokes system): 1.15166e-06, 8.87349e-17, 1.26557e-16, 1.43204e-16, 8.69635e-17, 1.43204e-16, 1.43204e-16, 1.39448e-16, 1.39812e-16, 1.39448e-16, 1.39448e-16, 0
Relative nonlinear residual (total system) after nonlinear iteration 1: 1.15166e-06


Postprocessing:
Compositions min/max/mass: 0.02515/0.1554/5.079e+09 // 0.00301/0.09499/9.701e+08 // 0.4/0.4/1.6e+10 // 0.3/0.3/1.2e+10 // 0.2/0.2/8e+09 // 0.1/0.1/4e+09 // 0.01006/0.06214/2.032e+09 // 0.007546/0.04661/1.524e+09 // 0.00503/0.03107/1.016e+09 // 0.002515/0.01554/5.079e+08

Termination requested by criterion: end time


+----------------------------------------------+------------+------------+
| Total wallclock time elapsed since start | 0.298s | |
| | | |
| Section | no. calls | wall time | % of total |
+----------------------------------+-----------+------------+------------+
| Assemble Stokes system rhs | 2 | 0.00562s | 1.9% |
| Assemble composition system | 20 | 0.071s | 24% |
| Assemble temperature system | 2 | 0.00398s | 1.3% |
| Build Stokes preconditioner | 2 | 0.000314s | 0.11% |
| Build composition preconditioner | 15 | 0.00126s | 0.42% |
| Build temperature preconditioner | 2 | 0.000281s | 0% |
| Initialization | 1 | 0.122s | 41% |
| Postprocessing | 2 | 0.00138s | 0.46% |
| Setup dof systems | 1 | 0.00306s | 1% |
| Setup initial conditions | 1 | 0.00475s | 1.6% |
| Setup matrices | 1 | 0.0031s | 1% |
| Solve Stokes system | 2 | 0.00204s | 0.69% |
| Solve composition reactions | 1 | 0.0676s | 23% |
| Solve composition system | 15 | 0.00337s | 1.1% |
| Solve temperature system | 2 | 0.000535s | 0.18% |
+----------------------------------+-----------+------------+------------+

-- Total wallclock time elapsed including restarts: 0s
-----------------------------------------------------------------------------
-- For information on how to cite ASPECT, see:
-- https://aspect.geodynamics.org/citing.html?ver=3.1.0-pre&mf=1&sha=986c9a918&src=code
-----------------------------------------------------------------------------
Loading