Skip to content

Commit ecfb751

Browse files
gassmoellergassmoeller
authored
Merge pull request #6161 from danieldouglas92/refactor_reactive_fluid_transport_model
Refactor the Tian Approximation in the Reactive Fluid Transport Model
2 parents 735206e + 215f207 commit ecfb751

File tree

11 files changed

+471
-209
lines changed

11 files changed

+471
-209
lines changed

Diff for: ‎contrib/utilities/update_scripts/prm_files/move_tian_reaction_model.py

+66
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,66 @@
1+
#!/usr/bin/env python3
2+
# -*- coding: utf-8 -*-
3+
4+
""" This script reformats the given .prm files to move the Tian 2019 reaction
5+
mode parameters to the correct subsection.
6+
"""
7+
8+
import sys
9+
import os
10+
import re
11+
import argparse
12+
13+
__author__ = 'The authors of the ASPECT code'
14+
__copyright__ = 'Copyright 2024, ASPECT'
15+
__license__ = 'GNU GPL 2 or later'
16+
17+
# Add the ASPECT root directory to the path so we can import from the aspect_data module
18+
sys.path.append(os.path.join(os.path.dirname(__file__), '..', '..', '..'))
19+
import python.scripts.aspect_input as aspect
20+
21+
22+
23+
def create_tian2019_reaction_subsection(parameters):
24+
""" Move the Tian 2019 reaction parameters to their own subsection. """
25+
26+
27+
parameters_to_move = ["Maximum weight percent water in sediment", \
28+
"Maximum weight percent water in MORB", \
29+
"Maximum weight percent water in gabbro", \
30+
"Maximum weight percent water in peridotite"]
31+
32+
# Collect existing parameters and delete old entries
33+
reactive_fluid_params = dict({})
34+
if "Reactive Fluid Transport Model" in parameters["Material model"]["value"]:
35+
reactive_fluid_params = parameters["Material model"]["value"]["Reactive Fluid Transport Model"]
36+
for param in parameters_to_move:
37+
if "Tian 2019 model" not in reactive_fluid_params["value"]:
38+
reactive_fluid_params["value"]["Tian 2019 model"] = {"comment": "", "value" : dict({}), "type": "subsection"}
39+
40+
if param in reactive_fluid_params["value"]:
41+
reactive_fluid_params["value"]["Tian 2019 model"]["value"][param] = reactive_fluid_params["value"][param]
42+
del reactive_fluid_params["value"][param]
43+
44+
return parameters
45+
46+
47+
48+
def main(input_file, output_file):
49+
"""Echo the input arguments to standard output"""
50+
parameters = aspect.read_parameter_file(input_file)
51+
52+
parameters = create_tian2019_reaction_subsection(parameters)
53+
54+
aspect.write_parameter_file(parameters, output_file)
55+
56+
57+
58+
if __name__ == '__main__':
59+
parser = argparse.ArgumentParser(
60+
prog='ASPECT .prm file reformatter',
61+
description='Reformats ASPECT .prm files to follow our general formatting guidelines. See the documentation of this script for details.')
62+
parser.add_argument('input_file', type=str, help='The .prm file to reformat')
63+
parser.add_argument('output_file', type=str, help='The .prm file to write the reformatted file to')
64+
args = parser.parse_args()
65+
66+
sys.exit(main(args.input_file, args.output_file))

Diff for: ‎cookbooks/tian_parameterization_kinematic_slab/coupled-two-phase-tian-parameterization-kinematic-slab.prm

+6-4
Original file line numberDiff line numberDiff line change
@@ -229,10 +229,12 @@ subsection Material model
229229
# values to encourage water to hydrate the overlying mantle. The polynomials defined
230230
# in Tian et al., 2019 also reach very large values at low P-T conditions, and so limiting
231231
# the weight percent to reasonable values is recommended.
232-
set Maximum weight percent water in peridotite = 2
233-
set Maximum weight percent water in gabbro = 1
234-
set Maximum weight percent water in MORB = 2
235-
set Maximum weight percent water in sediment = 3
232+
subsection Tian 2019 model
233+
set Maximum weight percent water in peridotite = 2
234+
set Maximum weight percent water in gabbro = 1
235+
set Maximum weight percent water in MORB = 2
236+
set Maximum weight percent water in sediment = 3
237+
end
236238
end
237239

238240
subsection Visco Plastic

Diff for: ‎doc/modules/changes/20250206_danieldouglas92

+5
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
Changed: Separated out the tian2019 solubility
2+
reaction model into it's own module independent
3+
of the reactive fluid transport material model.
4+
<br>
5+
(Daniel Douglas, 2025/02/06)

Diff for: ‎include/aspect/material_model/reaction_model/tian2019_solubility.h

+155
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
/*
2+
Copyright (C) 2024 by the authors of the ASPECT code.
3+
4+
This file is part of ASPECT.
5+
6+
ASPECT is free software; you can redistribute it and/or modify
7+
it under the terms of the GNU General Public License as published by
8+
the Free Software Foundation; either version 2, or (at your option)
9+
any later version.
10+
11+
ASPECT is distributed in the hope that it will be useful,
12+
but WITHOUT ANY WARRANTY; without even the implied warranty of
13+
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14+
GNU General Public License for more details.
15+
16+
You should have received a copy of the GNU General Public License
17+
along with ASPECT; see the file LICENSE. If not see
18+
<http://www.gnu.org/licenses/>.
19+
*/
20+
21+
#ifndef _aspect_material_reaction_model_melt_tian2019_solubility_h
22+
#define _aspect_material_reaction_model_melt_tian2019_solubility_h
23+
24+
#include <aspect/material_model/interface.h>
25+
#include <aspect/simulator_access.h>
26+
#include <aspect/melt.h>
27+
#include <aspect/utilities.h>
28+
29+
30+
namespace aspect
31+
{
32+
namespace MaterialModel
33+
{
34+
using namespace dealii;
35+
36+
namespace ReactionModel
37+
{
38+
39+
/**
40+
* A melt model that calculates the solubility of water according to
41+
* parameterized phase diagrams for four lithologies:
42+
* 1) sediment
43+
* 2) mid-ocean ridge basalt (MORB)
44+
* 3) gabbro
45+
* 4) peridotite
46+
* from Tian, 2019 https://doi.org/10.1029/2019GC008488.
47+
*
48+
* These functions can be used in the calculation of reactive fluid transport
49+
* of water.
50+
*
51+
* @ingroup ReactionModel
52+
*/
53+
template <int dim>
54+
class Tian2019Solubility : public ::aspect::SimulatorAccess<dim>
55+
{
56+
public:
57+
58+
/**
59+
* Compute the free fluid fraction that is present in the material based on the
60+
* fluid content of the material and the fluid solubility for the given input conditions.
61+
* @p in and @p melt_fraction need to have the same size.
62+
*
63+
* @param in Object that contains the current conditions.
64+
* @param porosity_idx the index of the "porosity" composition
65+
* @param q the quadrature point index
66+
*/
67+
double
68+
melt_fraction (const MaterialModel::MaterialModelInputs<dim> &in,
69+
const unsigned int porosity_idx,
70+
unsigned int q) const;
71+
72+
/**
73+
* Compute the maximum allowed bound water content at the input
74+
* pressure and temperature conditions. This is used to determine
75+
* how free water interacts with the solid phase.
76+
* @param in Object that contains the current conditions.
77+
* @param q the quadrature point index
78+
*/
79+
std::vector<double> tian_equilibrium_bound_water_content(const MaterialModel::MaterialModelInputs<dim> &in,
80+
unsigned int q) const;
81+
82+
/**
83+
* Declare the parameters this function takes through input files.
84+
*/
85+
static
86+
void
87+
declare_parameters (ParameterHandler &prm);
88+
89+
/**
90+
* Read the parameters from the parameter file.
91+
*/
92+
void
93+
parse_parameters (ParameterHandler &prm);
94+
95+
private:
96+
97+
/**
98+
* The maximum water content for each of the 4 rock types in the tian approximation
99+
* method. These are important for keeping the polynomial bounded within reasonable
100+
* values.
101+
*/
102+
double tian_max_peridotite_water;
103+
double tian_max_gabbro_water;
104+
double tian_max_MORB_water;
105+
double tian_max_sediment_water;
106+
107+
/**
108+
*
109+
* The following coefficients are taken from a publication from Tian et al., 2019, and can be found
110+
* in Table 3 (Gabbro), Table B1 (MORB), Table B2 (Sediments) and Table B3 (peridotite).
111+
* LR refers to the effective enthalpy change for devolatilization reactions,
112+
* csat is the saturated mass fraction of water in the solid, and Td is the
113+
* onset temperature of devolatilization for water.
114+
*/
115+
std::vector<double> LR_peridotite_poly_coeffs {-19.0609, 168.983, -630.032, 1281.84, -1543.14, 1111.88, -459.142, 95.4143, 1.97246};
116+
std::vector<double> csat_peridotite_poly_coeffs {0.00115628, 2.42179};
117+
std::vector<double> Td_peridotite_poly_coeffs {-15.4627, 94.9716, 636.603};
118+
119+
std::vector<double> LR_gabbro_poly_coeffs {-1.81745, 7.67198, -10.8507, 5.09329, 8.14519};
120+
std::vector<double> csat_gabbro_poly_coeffs {-0.0176673, 0.0893044, 1.52732};
121+
std::vector<double> Td_gabbro_poly_coeffs {-1.72277, 20.5898, 637.517};
122+
123+
std::vector<double> LR_MORB_poly_coeffs {-1.78177, 7.50871, -10.4840, 5.19725, 7.96365};
124+
std::vector<double> csat_MORB_poly_coeffs {0.0102725, -0.115390, 0.324452, 1.41588};
125+
std::vector<double> Td_MORB_poly_coeffs {-3.81280, 22.7809, 638.049};
126+
127+
std::vector<double> LR_sediment_poly_coeffs {-2.03283, 10.8186, -21.2119, 18.3351, -6.48711, 8.32459};
128+
std::vector<double> csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867};
129+
std::vector<double> Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898};
130+
131+
/**
132+
* The polynomials breakdown above certain pressures, 10 GPa for peridotite, 26 GPa for gabbro, 16 GPa for MORB,
133+
* and 50 GPa for sediment. These cutoff pressures were determined by extending the pressure range in Tian et al. (2019)
134+
* and observing where the maximum allowed water contents jump towards infinite values.
135+
*/
136+
const std::array<double,4 > pressure_cutoffs {{10, 26, 16, 50}};
137+
138+
std::vector<std::vector<double>> devolatilization_enthalpy_changes {LR_peridotite_poly_coeffs, LR_gabbro_poly_coeffs, \
139+
LR_MORB_poly_coeffs, LR_sediment_poly_coeffs
140+
};
141+
142+
std::vector<std::vector<double>> water_mass_fractions {csat_peridotite_poly_coeffs, csat_gabbro_poly_coeffs, \
143+
csat_MORB_poly_coeffs, csat_sediment_poly_coeffs
144+
};
145+
146+
std::vector<std::vector<double>> devolatilization_onset_temperatures {Td_peridotite_poly_coeffs, Td_gabbro_poly_coeffs, \
147+
Td_MORB_poly_coeffs, Td_sediment_poly_coeffs
148+
};
149+
};
150+
}
151+
152+
}
153+
}
154+
155+
#endif

Diff for: ‎include/aspect/material_model/reactive_fluid_transport.h

+6-69
Original file line numberDiff line numberDiff line change
@@ -28,11 +28,8 @@
2828
#include <aspect/geometry_model/interface.h>
2929

3030
#include <aspect/melt.h>
31-
#include <aspect/utilities.h>
32-
#include <aspect/geometry_model/interface.h>
3331
#include <aspect/material_model/reaction_model/katz2003_mantle_melting.h>
34-
35-
32+
#include <aspect/material_model/reaction_model/tian2019_solubility.h>
3633

3734
namespace aspect
3835
{
@@ -63,18 +60,6 @@ namespace aspect
6360
* @}
6461
*/
6562

66-
67-
/**
68-
* Compute the maximum allowed bound water content at the input
69-
* pressure and temperature conditions. This is used to determine
70-
* how free water interacts with the solid phase.
71-
* @param in Object that contains the current conditions.
72-
* @param q unsigned int from 0-3 indexing which rock phase the equilbrium
73-
* bound water content is being calculated for
74-
*/
75-
std::vector<double> tian_equilibrium_bound_water_content(const MaterialModel::MaterialModelInputs<dim> &in,
76-
unsigned int q) const;
77-
7863
/**
7964
* Compute the free fluid fraction that can be present in the material based on the
8065
* fluid content of the material and the fluid solubility for the given input conditions.
@@ -159,64 +144,16 @@ namespace aspect
159144
*/
160145
double fluid_reaction_time_scale;
161146

162-
/**
163-
* The maximum water content for each of the 4 rock types in the tian approximation
164-
* method. These are important for keeping the polynomial bounded within reasonable
165-
* values.
166-
*/
167-
double tian_max_peridotite_water;
168-
double tian_max_gabbro_water;
169-
double tian_max_MORB_water;
170-
double tian_max_sediment_water;
171-
172-
/**
173-
*
174-
* The following coefficients are taken from a publication from Tian et al., 2019, and can be found
175-
* in Table 3 (Gabbro), Table B1 (MORB), Table B2 (Sediments) and Table B3 (peridotite).
176-
* LR refers to the effective enthalpy change for devolatilization reactions,
177-
* csat is the saturated mass fraction of water in the solid, and Td is the
178-
* onset temperature of devolatilization for water.
179-
*/
180-
std::vector<double> LR_peridotite_poly_coeffs {-19.0609, 168.983, -630.032, 1281.84, -1543.14, 1111.88, -459.142, 95.4143, 1.97246};
181-
std::vector<double> csat_peridotite_poly_coeffs {0.00115628, 2.42179};
182-
std::vector<double> Td_peridotite_poly_coeffs {-15.4627, 94.9716, 636.603};
183-
184-
std::vector<double> LR_gabbro_poly_coeffs {-1.81745, 7.67198, -10.8507, 5.09329, 8.14519};
185-
std::vector<double> csat_gabbro_poly_coeffs {-0.0176673, 0.0893044, 1.52732};
186-
std::vector<double> Td_gabbro_poly_coeffs {-1.72277, 20.5898, 637.517};
187-
188-
std::vector<double> LR_MORB_poly_coeffs {-1.78177, 7.50871, -10.4840, 5.19725, 7.96365};
189-
std::vector<double> csat_MORB_poly_coeffs {0.0102725, -0.115390, 0.324452, 1.41588};
190-
std::vector<double> Td_MORB_poly_coeffs {-3.81280, 22.7809, 638.049};
191-
192-
std::vector<double> LR_sediment_poly_coeffs {-2.03283, 10.8186, -21.2119, 18.3351, -6.48711, 8.32459};
193-
std::vector<double> csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867};
194-
std::vector<double> Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898};
195-
196-
/**
197-
* The polynomials breakdown above certain pressures, 10 GPa for peridotite, 26 GPa for gabbro, 16 GPa for MORB,
198-
* and 50 GPa for sediment. These cutoff pressures were determined by extending the pressure range in Tian et al. (2019)
199-
* and observing where the maximum allowed water contents jump towards infinite values.
200-
*/
201-
const std::vector<double> pressure_cutoffs {10, 26, 16, 50};
202-
203-
std::vector<std::vector<double>> devolatilization_enthalpy_changes {LR_peridotite_poly_coeffs, LR_gabbro_poly_coeffs, \
204-
LR_MORB_poly_coeffs, LR_sediment_poly_coeffs
205-
};
206-
207-
std::vector<std::vector<double>> water_mass_fractions {csat_peridotite_poly_coeffs, csat_gabbro_poly_coeffs, \
208-
csat_MORB_poly_coeffs, csat_sediment_poly_coeffs
209-
};
210-
211-
std::vector<std::vector<double>> devolatilization_onset_temperatures {Td_peridotite_poly_coeffs, Td_gabbro_poly_coeffs, \
212-
Td_MORB_poly_coeffs, Td_sediment_poly_coeffs
213-
};
214-
215147
/*
216148
* Object for computing Katz 2003 melt parameters
217149
*/
218150
ReactionModel::Katz2003MantleMelting<dim> katz2003_model;
219151

152+
/*
153+
* Object for computing Tian 2019 parameterized solubility parameters
154+
*/
155+
ReactionModel::Tian2019Solubility<dim> tian2019_model;
156+
220157
/**
221158
* Enumeration for selecting which type of scheme to use for
222159
* reactions between fluids and solids. The available

0 commit comments

Comments
 (0)