Skip to content

Commit 0b13719

Browse files
danieldouglas92danieldouglas92
committed
refactor the tian approximation in the reactive fluid transport model
1 parent 8b5abd1 commit 0b13719

File tree

9 files changed

+400
-209
lines changed

9 files changed

+400
-209
lines changed

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

Lines changed: 6 additions & 4 deletions
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

‎include/aspect/material_model/reaction_model/tian2019_solubility.h

Lines changed: 155 additions & 0 deletions
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

‎include/aspect/material_model/reactive_fluid_transport.h

Lines changed: 6 additions & 69 deletions
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
{
@@ -64,18 +61,6 @@ namespace aspect
6461
* @}
6562
*/
6663

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

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

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

0 commit comments

Comments
 (0)