diff --git a/include/aspect/material_model/reaction_model/tian2019_solubility.h b/include/aspect/material_model/reaction_model/tian2019_solubility.h index e27b4e68802..3957960a9c9 100644 --- a/include/aspect/material_model/reaction_model/tian2019_solubility.h +++ b/include/aspect/material_model/reaction_model/tian2019_solubility.h @@ -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: /** diff --git a/source/material_model/reaction_model/tian2019_solubility.cc b/source/material_model/reaction_model/tian2019_solubility.cc index 640fe272856..d7770bc4b63 100644 --- a/source/material_model/reaction_model/tian2019_solubility.cc +++ b/source/material_model/reaction_model/tian2019_solubility.cc @@ -147,6 +147,13 @@ namespace aspect void Tian2019Solubility::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."); @@ -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"); diff --git a/source/material_model/reactive_fluid_transport.cc b/source/material_model/reactive_fluid_transport.cc index c719593690c..76be2890c20 100644 --- a/source/material_model/reactive_fluid_transport.cc +++ b/source/material_model/reactive_fluid_transport.cc @@ -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)) @@ -252,6 +252,32 @@ namespace aspect std::vector 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; qintrospection().compositional_index_for_name("bound_fluid"); @@ -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; @@ -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 / @@ -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; + } } } } diff --git a/tests/tian_track_fluid_sources.prm b/tests/tian_track_fluid_sources.prm new file mode 100644 index 00000000000..c6cf0663e27 --- /dev/null +++ b/tests/tian_track_fluid_sources.prm @@ -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 diff --git a/tests/tian_track_fluid_sources/screen-output b/tests/tian_track_fluid_sources/screen-output new file mode 100644 index 00000000000..baac228113b --- /dev/null +++ b/tests/tian_track_fluid_sources/screen-output @@ -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 +----------------------------------------------------------------------------- diff --git a/tests/tian_track_fluid_sources/statistics b/tests/tian_track_fluid_sources/statistics new file mode 100644 index 00000000000..51a3b1de10e --- /dev/null +++ b/tests/tian_track_fluid_sources/statistics @@ -0,0 +1,54 @@ +# 1: Time step number +# 2: Time (years) +# 3: Time step size (years) +# 4: Number of mesh cells +# 5: Number of Stokes degrees of freedom +# 6: Number of temperature degrees of freedom +# 7: Number of degrees of freedom for all compositions +# 8: Number of nonlinear iterations +# 9: Iterations for temperature solver +# 10: Iterations for composition solver 1 +# 11: Iterations for composition solver 2 +# 12: Iterations for composition solver 3 +# 13: Iterations for composition solver 4 +# 14: Iterations for composition solver 5 +# 15: Iterations for composition solver 6 +# 16: Iterations for composition solver 7 +# 17: Iterations for composition solver 8 +# 18: Iterations for composition solver 9 +# 19: Iterations for composition solver 10 +# 20: Iterations for Stokes solver +# 21: Velocity iterations in Stokes preconditioner +# 22: Schur complement iterations in Stokes preconditioner +# 23: Minimal value for composition porosity +# 24: Maximal value for composition porosity +# 25: Global mass for composition porosity +# 26: Minimal value for composition bound_fluid +# 27: Maximal value for composition bound_fluid +# 28: Global mass for composition bound_fluid +# 29: Minimal value for composition peridotite +# 30: Maximal value for composition peridotite +# 31: Global mass for composition peridotite +# 32: Minimal value for composition gabbro +# 33: Maximal value for composition gabbro +# 34: Global mass for composition gabbro +# 35: Minimal value for composition MORB +# 36: Maximal value for composition MORB +# 37: Global mass for composition MORB +# 38: Minimal value for composition sediment +# 39: Maximal value for composition sediment +# 40: Global mass for composition sediment +# 41: Minimal value for composition peridotite_porosity +# 42: Maximal value for composition peridotite_porosity +# 43: Global mass for composition peridotite_porosity +# 44: Minimal value for composition gabbro_porosity +# 45: Maximal value for composition gabbro_porosity +# 46: Global mass for composition gabbro_porosity +# 47: Minimal value for composition MORB_porosity +# 48: Maximal value for composition MORB_porosity +# 49: Global mass for composition MORB_porosity +# 50: Minimal value for composition sediment_porosity +# 51: Maximal value for composition sediment_porosity +# 52: Global mass for composition sediment_porosity +0 0.000000000000e+00 0.000000000000e+00 16 187 81 810 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.00000000e+00 0.00000000e+00 0.00000000e+00 1.10000000e-01 1.10000000e-01 4.40000000e+09 4.00000000e-01 4.00000000e-01 1.60000000e+10 3.00000000e-01 3.00000000e-01 1.20000000e+10 2.00000000e-01 2.00000000e-01 8.00000000e+09 1.00000000e-01 1.00000000e-01 4.00000000e+09 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 0.00000000e+00 +1 1.000000000000e+03 1.000000000000e+03 16 187 81 810 1 6 0 0 0 0 0 0 0 0 0 0 0 0 0 2.51518787e-02 1.55352795e-01 5.07898344e+09 3.01021900e-03 9.49916812e-02 9.70052691e+08 4.00000000e-01 4.00000000e-01 1.60000000e+10 3.00000000e-01 3.00000000e-01 1.20000000e+10 2.00000000e-01 2.00000000e-01 8.00000000e+09 1.00000000e-01 1.00000000e-01 4.00000000e+09 1.00607515e-02 6.21411181e-02 2.03159337e+09 7.54556360e-03 4.66058386e-02 1.52369503e+09 5.03037573e-03 3.10705590e-02 1.01579669e+09 2.51518787e-03 1.55352795e-02 5.07898344e+08