Skip to content

Commit

Permalink
Address comments
Browse files Browse the repository at this point in the history
  • Loading branch information
@danieldouglas92
danieldouglas92 committed Feb 25, 2024
1 parent e0ac649 commit d04bdca
Show file tree
Hide file tree
Showing 15 changed files with 137 additions and 92 deletions.
5 changes: 3 additions & 2 deletions doc/modules/changes/20230714_danieldouglas92 → doc/modules/changes/20240224_danieldouglas92
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
Implemented a new fluid-rock interaction scheme for the reactive fluid
transport material model based on a parametrization by Tian et al., 2018.
transport material model based on a parametrization by Tian et al., 2019 G3,
https://doi.org/10.1029/2019GC008488.
<br>
(Daniel Douglas, 2023/07/14)
(Daniel Douglas, 2024/02/24)
8 changes: 4 additions & 4 deletions include/aspect/material_model/reactive_fluid_transport.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,10 +141,10 @@ namespace aspect

// Time scale for fluid release and absorption.
double fluid_reaction_time_scale;
double max_peridotite_water;
double max_gabbro_water;
double max_MORB_water;
double max_sediment_water;
double tian_max_peridotite_water;
double tian_max_gabbro_water;
double tian_max_MORB_water;
double tian_max_sediment_water;

/**
* Enumeration for selecting which type of scheme to use for
Expand Down
65 changes: 38 additions & 27 deletions source/material_model/reactive_fluid_transport.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,9 +84,11 @@ namespace aspect
std::vector<std::vector<double>> LR_all_poly_coeffs {LR_peridotite_poly_coeffs, LR_gabbro_poly_coeffs, \
LR_MORB_poly_coeffs, LR_sediment_poly_coeffs
};

std::vector<std::vector<double>> csat_all_poly_coeffs {csat_peridotite_poly_coeffs, csat_gabbro_poly_coeffs, \
csat_MORB_poly_coeffs, csat_sediment_poly_coeffs
};

std::vector<std::vector<double>> Td_all_poly_coeffs {Td_peridotite_poly_coeffs, Td_gabbro_poly_coeffs, \
Td_MORB_poly_coeffs, Td_sediment_poly_coeffs
};
Expand All @@ -96,8 +98,8 @@ namespace aspect
std::vector<double> csat_values {0, 0, 0, 0};
std::vector<double> Td_values {0, 0, 0, 0};

// Loop over the four rock types i (peridotite, gabbro, MORB, sediment) and the polynomial
// coefficients j to fill the vectors defined above. The polynomials for LR are defined in
// Loop over the four rock types (peridotite, gabbro, MORB, sediment) and the polynomial
// coefficients to fill the vectors defined above. The polynomials for LR are defined in
// equations 13, B2, B10, and B18. csat polynomials are defined in equations 14, B1, B9, and B17.
// Td polynomials are defined in equations 15, B3, B11, and B19.
for (unsigned int i = 0; i<LR_all_poly_coeffs.size(); ++i)
Expand All @@ -123,7 +125,7 @@ namespace aspect
std::vector<double> eq_bound_water_content;

// Define the maximum bound water content allowed for the four different rock compositions
std::vector<double> max_bound_water_content = {max_peridotite_water, max_gabbro_water, max_MORB_water, max_sediment_water};
std::vector<double> max_bound_water_content = {tian_max_peridotite_water, tian_max_gabbro_water, tian_max_MORB_water, tian_max_sediment_water};

// Loop over all rock compositions and fill the equilibrium bound water content, divide by 100 to convert
// from percentage to fraction (equation 1)
Expand Down Expand Up @@ -152,7 +154,6 @@ namespace aspect
{
case no_reaction:
{

// No reactions occur between the solid and fluid phases,
// and the fluid volume fraction (stored in the melt_fractions
// vector) is equal to the porosity.
Expand All @@ -161,34 +162,38 @@ namespace aspect
}
case zero_solubility:
{
// A fluid-rock reaction model where no reactions occur.
// The melt (fluid) fraction at any point is equal
// to the sum of the bound fluid content and porosity,
// with the latter determined by the assigned initial
// porosity, fluid boundary conditions, and fluid
// transport through the model.
// The fluid volume fraction in equilibrium with the solid
// at any point (stored in the melt_fractions vector) is
// equal to the sum of the bound fluid content and porosity.
melt_fractions[q] = in.composition[q][bound_fluid_idx] + in.composition[q][porosity_idx];
break;
}
case tian_approximation:
{
// The bound fluid content is calculated using parametrized phase
// diagrams for four different rock types: sediment, MORB, gabbro, and
// peridotite.
const unsigned int sediment_idx = this->introspection().compositional_index_for_name("sediment");
const unsigned int MORB_idx = this->introspection().compositional_index_for_name("MORB");
const unsigned int gabbro_idx = this->introspection().compositional_index_for_name("gabbro");
const unsigned int peridotite_idx = this->introspection().compositional_index_for_name("peridotite");

// Initialize a vector that stores the compositions which tracks the four different rock compositions,
// and these compositions are tracked as mass fractions
std::vector<double> tracked_rock_compositions;
tracked_rock_compositions.push_back(in.composition[q][peridotite_idx]);
tracked_rock_compositions.push_back(in.composition[q][gabbro_idx]);
tracked_rock_compositions.push_back(in.composition[q][MORB_idx]);
tracked_rock_compositions.push_back(in.composition[q][sediment_idx]);
// Initialize a vector that stores the compositions (mass fractions) for
// the four different rock compositions,
std::vector<double> tracked_rock_mass_fractions;
tracked_rock_mass_fractions.push_back(in.composition[q][peridotite_idx]);
tracked_rock_mass_fractions.push_back(in.composition[q][gabbro_idx]);
tracked_rock_mass_fractions.push_back(in.composition[q][MORB_idx]);
tracked_rock_mass_fractions.push_back(in.composition[q][sediment_idx]);

// The bound water content (water within the solid phase) for the four different rock types
std::vector<double> tian_eq_bound_water_content = tian_equilibrium_bound_water_content(in, q);

double average_eq_bound_water_content = MaterialUtilities::average_value (tracked_rock_compositions, tian_eq_bound_water_content, MaterialUtilities::arithmetic);
// average the water content between the four different rock types
double average_eq_bound_water_content = MaterialUtilities::average_value (tracked_rock_mass_fractions, tian_eq_bound_water_content, MaterialUtilities::arithmetic);

// The fluid volume fraction (stored in the melt_fractions vector) is equal to the porosity and the change
// in bound fluid content (current bound fluid - updated average bound fluid).
melt_fractions[q] = std::max(in.composition[q][bound_fluid_idx] + in.composition[q][porosity_idx] - average_eq_bound_water_content, 0.0);
break;
}
Expand Down Expand Up @@ -227,9 +232,10 @@ namespace aspect
typename Interface<dim>::MaterialModelOutputs &out) const
{
base_model->evaluate(in,out);
// Modify the viscosity from the base model based on the presence of fluid.

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))
{
// Scale the base model viscosity value based on the porosity.
Expand Down Expand Up @@ -356,16 +362,16 @@ namespace aspect
"computed. If the model does not use operator splitting, this parameter is not used. "
"Units: yr or s, depending on the ``Use years "
"in output instead of seconds'' parameter.");
prm.declare_entry ("Maximum weight percent sediment", "3",
prm.declare_entry ("Maximum weight percent water in sediment", "3",
Patterns::Double (0),
"The maximum allowed weight percent that the sediment composition can hold.");
prm.declare_entry ("Maximum weight percent MORB", "2",
prm.declare_entry ("Maximum weight percent water in MORB", "2",
Patterns::Double (0),
"The maximum allowed weight percent that the sediment composition can hold.");
prm.declare_entry ("Maximum weight percent gabbro", "1",
prm.declare_entry ("Maximum weight percent water in gabbro", "1",
Patterns::Double (0),
"The maximum allowed weight percent that the sediment composition can hold.");
prm.declare_entry ("Maximum weight percent peridotite", "8",
prm.declare_entry ("Maximum weight percent water in peridotite", "8",
Patterns::Double (0),
"The maximum allowed weight percent that the sediment composition can hold.");
prm.declare_entry ("Fluid-solid reaction scheme", "no reaction",
Expand Down Expand Up @@ -405,10 +411,10 @@ namespace aspect
fluid_compressibility = prm.get_double ("Fluid compressibility");
fluid_reaction_time_scale = prm.get_double ("Fluid reaction time scale for operator splitting");

max_peridotite_water = prm.get_double ("Maximum weight percent peridotite");
max_gabbro_water = prm.get_double ("Maximum weight percent gabbro");
max_MORB_water = prm.get_double ("Maximum weight percent MORB");
max_sediment_water = prm.get_double ("Maximum weight percent sediment");
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");
tian_max_sediment_water = prm.get_double ("Maximum weight percent water in sediment");

// Create the base model and initialize its SimulatorAccess base
// class; it will get a chance to read its parameters below after we
Expand Down Expand Up @@ -493,6 +499,11 @@ namespace aspect
// After parsing the parameters for this model, parse parameters related to the base model.
base_model->parse_parameters(prm);
this->model_dependence = base_model->get_model_dependence();
if (fluid_solid_reaction_scheme == zero_solubility)
{
AssertThrow(this->get_material_model().is_compressible() == false,
ExcMessage("The Fluid-reaction scheme zero solubility must be used with an incompressible base model."));
}
}


Expand Down
6 changes: 3 additions & 3 deletions tests/tian_MORB.prm
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# 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 reproduces that phase diagram, with an important
# wt% water that can be stored in a MORB as shown in Figure 5 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.
Expand All @@ -20,6 +20,6 @@ end

subsection Material model
subsection Reactive Fluid Transport Model
set Maximum weight percent MORB = 5.3
set Maximum weight percent water in MORB = 5.3
end
end
15 changes: 13 additions & 2 deletions tests/tian_MORB/screen-output
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
Number of active cells: 16 (on 3 levels)
Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)

Expand All @@ -17,7 +23,7 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)


Postprocessing:
Writing graphical output: output/solution/solution-00000
Writing graphical output: output-tian_MORB/solution/solution-00000
Compositions min/max/mass: 0/0/0 // 0.11/0.11/4.4e+09 // 0/0/0 // 0/0/0 // 1/1/4e+10 // 0/0/0
RMS, max velocity: 5.54e-21 m/year, 1.26e-20 m/year

Expand All @@ -38,11 +44,16 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)


Postprocessing:
Writing graphical output: output/solution/solution-00001
Writing graphical output: output-tian_MORB/solution/solution-00001
Compositions min/max/mass: 0.057/0.1095/3.858e+09 // 0.0005023/0.053/5.422e+08 // 0/0/0 // 0/0/0 // 1/1/4e+10 // 0/0/0
RMS, max velocity: 5.54e-21 m/year, 1.26e-20 m/year

Termination requested by criterion: end time


+----------------------------------------------+------------+------------+
+----------------------------------+-----------+------------+------------+
+----------------------------------+-----------+------------+------------+

-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
4 changes: 2 additions & 2 deletions tests/tian_MORB/statistics
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,5 +37,5 @@
# 37: Global mass for composition sediment
# 38: RMS velocity (m/year)
# 39: Max. velocity (m/year)
0 0.000000000000e+00 0.000000000000e+00 16 187 81 486 1 0 0 0 0 0 0 0 12 14 90 output/solution/solution-00000 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.10000000e-01 1.10000000e-01 4.40000000e+09 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 4.00000000e+10 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.54090700e-21 1.25942419e-20
1 1.000000000000e+03 1.000000000000e+03 16 187 81 486 1 0 0 0 0 0 0 0 4294967295 0 0 output/solution/solution-00001 5.70000000e-02 1.09497728e-01 3.85780540e+09 5.02272469e-04 5.30000000e-02 5.42194599e+08 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 4.00000000e+10 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.54090700e-21 1.25942419e-20
0 0.000000000000e+00 0.000000000000e+00 16 187 81 486 1 0 0 0 0 0 0 0 12 14 90 output-tian_MORB/solution/solution-00000 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.10000000e-01 1.10000000e-01 4.40000000e+09 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 4.00000000e+10 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.54090700e-21 1.25942419e-20
1 1.000000000000e+03 1.000000000000e+03 16 187 81 486 1 0 0 0 0 0 0 0 4294967295 0 0 output-tian_MORB/solution/solution-00001 5.70000000e-02 1.09497728e-01 3.85780540e+09 5.02272469e-04 5.30000000e-02 5.42194599e+08 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.00000000e+00 1.00000000e+00 4.00000000e+10 0.00000000e+00 0.00000000e+00 0.00000000e+00 5.54090700e-21 1.25942419e-20
6 changes: 3 additions & 3 deletions tests/tian_gabbro.prm
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# 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 reproduces that phase diagram, with an important
# wt% water that can be stored in a gabbro as shown in Figure 2 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.
Expand All @@ -20,6 +20,6 @@ end

subsection Material model
subsection Reactive Fluid Transport Model
set Maximum weight percent gabbro = 5.1
set Maximum weight percent water in gabbro = 5.1
end
end
15 changes: 13 additions & 2 deletions tests/tian_gabbro/screen-output
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
-----------------------------------------------------------------------------
Number of active cells: 16 (on 3 levels)
Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)

Expand All @@ -17,7 +23,7 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)


Postprocessing:
Writing graphical output: output/solution/solution-00000
Writing graphical output: output-tian_gabbro/solution/solution-00000
Compositions min/max/mass: 0/0/0 // 0.11/0.11/4.4e+09 // 0/0/0 // 1/1/4e+10 // 0/0/0 // 0/0/0
RMS, max velocity: 5.54e-21 m/year, 1.26e-20 m/year

Expand All @@ -38,11 +44,16 @@ Number of degrees of freedom: 989 (162+73+162+25+81+81+81+81+81+81+81)


Postprocessing:
Writing graphical output: output/solution/solution-00001
Writing graphical output: output-tian_gabbro/solution/solution-00001
Compositions min/max/mass: 0.059/0.1096/3.792e+09 // 0.0003505/0.051/6.078e+08 // 0/0/0 // 1/1/4e+10 // 0/0/0 // 0/0/0
RMS, max velocity: 5.54e-21 m/year, 1.26e-20 m/year

Termination requested by criterion: end time


+----------------------------------------------+------------+------------+
+----------------------------------+-----------+------------+------------+
+----------------------------------+-----------+------------+------------+

-----------------------------------------------------------------------------
-----------------------------------------------------------------------------

0 comments on commit d04bdca

Please sign in to comment.