Skip to content

Commit af3a3e0

Browse files
danieldouglas92danieldouglas92
committed
refactor the tian approximation in the reactive fluid transport model
1 parent b14f5fd commit af3a3e0

File tree

9 files changed

+397
-197
lines changed

9 files changed

+397
-197
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: 159 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,159 @@
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_melt_tian2019_solubility_h
22+
#define _aspect_material_reaction_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 melt fraction and entropy change
41+
* according to the melting model for dry peridotite of Katz, 2003.
42+
* This also includes a computation of the latent heat of melting (if the latent heat
43+
* heating model is active).
44+
*
45+
* These functions can be used in the calculation of melting and melt transport
46+
* in the melt_simple material model and can be extended to other material models
47+
*
48+
* @ingroup ReactionModel
49+
*/
50+
template <int dim>
51+
class Tian2019Solubility : public ::aspect::SimulatorAccess<dim>
52+
{
53+
public:
54+
// constructor
55+
Tian2019Solubility();
56+
57+
/**
58+
* Declare the parameters this function takes through input files.
59+
*/
60+
static
61+
void
62+
declare_parameters (ParameterHandler &prm);
63+
64+
/**
65+
* Read the parameters from the parameter file.
66+
*/
67+
void
68+
parse_parameters (ParameterHandler &prm);
69+
70+
71+
/**
72+
* Compute the free fluid fraction that can be present in the material based on the
73+
* fluid content of the material and the fluid solubility for the given input conditions.
74+
* @p in and @p melt_fractions need to have the same size.
75+
*
76+
* @param in Object that contains the current conditions.
77+
* @param melt_fractions Vector of doubles that is filled with the
78+
* allowable free fluid fraction for each given input conditions.
79+
*/
80+
double
81+
melt_fraction (const MaterialModel::MaterialModelInputs<dim> &in,
82+
const unsigned int porosity_idx,
83+
unsigned int q) const;
84+
85+
/**
86+
* Compute the maximum allowed bound water content at the input
87+
* pressure and temperature conditions. This is used to determine
88+
* how free water interacts with the solid phase.
89+
* @param in Object that contains the current conditions.
90+
* @param q unsigned int from 0-3 indexing which rock phase the equilbrium
91+
* bound water content is being calculated for
92+
*/
93+
std::vector<double> tian_equilibrium_bound_water_content(const MaterialModel::MaterialModelInputs<dim> &in,
94+
unsigned int q) const;
95+
96+
private:
97+
/**
98+
* Parameters for the solubility model of Tian et al., 2019.
99+
*/
100+
101+
/**
102+
* The maximum water content for each of the 4 rock types in the tian approximation
103+
* method. These are important for keeping the polynomial bounded within reasonable
104+
* values.
105+
*/
106+
double tian_max_peridotite_water;
107+
double tian_max_gabbro_water;
108+
double tian_max_MORB_water;
109+
double tian_max_sediment_water;
110+
111+
/**
112+
*
113+
* The following coefficients are taken from a publication from Tian et al., 2019, and can be found
114+
* in Table 3 (Gabbro), Table B1 (MORB), Table B2 (Sediments) and Table B3 (peridotite).
115+
* LR refers to the effective enthalpy change for devolatilization reactions,
116+
* csat is the saturated mass fraction of water in the solid, and Td is the
117+
* onset temperature of devolatilization for water.
118+
*/
119+
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};
120+
std::vector<double> csat_peridotite_poly_coeffs {0.00115628, 2.42179};
121+
std::vector<double> Td_peridotite_poly_coeffs {-15.4627, 94.9716, 636.603};
122+
123+
std::vector<double> LR_gabbro_poly_coeffs {-1.81745, 7.67198, -10.8507, 5.09329, 8.14519};
124+
std::vector<double> csat_gabbro_poly_coeffs {-0.0176673, 0.0893044, 1.52732};
125+
std::vector<double> Td_gabbro_poly_coeffs {-1.72277, 20.5898, 637.517};
126+
127+
std::vector<double> LR_MORB_poly_coeffs {-1.78177, 7.50871, -10.4840, 5.19725, 7.96365};
128+
std::vector<double> csat_MORB_poly_coeffs {0.0102725, -0.115390, 0.324452, 1.41588};
129+
std::vector<double> Td_MORB_poly_coeffs {-3.81280, 22.7809, 638.049};
130+
131+
std::vector<double> LR_sediment_poly_coeffs {-2.03283, 10.8186, -21.2119, 18.3351, -6.48711, 8.32459};
132+
std::vector<double> csat_sediment_poly_coeffs {-0.150662, 0.301807, 1.01867};
133+
std::vector<double> Td_sediment_poly_coeffs {2.83277, -24.7593, 85.9090, 524.898};
134+
135+
/**
136+
* The polynomials breakdown above certain pressures, 10 GPa for peridotite, 26 GPa for gabbro, 16 GPa for MORB,
137+
* and 50 GPa for sediment. These cutoff pressures were determined by extending the pressure range in Tian et al. (2019)
138+
* and observing where the maximum allowed water contents jump towards infinite values.
139+
*/
140+
const std::vector<double> pressure_cutoffs {10, 26, 16, 50};
141+
142+
std::vector<std::vector<double>> devolatilization_enthalpy_changes {LR_peridotite_poly_coeffs, LR_gabbro_poly_coeffs, \
143+
LR_MORB_poly_coeffs, LR_sediment_poly_coeffs
144+
};
145+
146+
std::vector<std::vector<double>> water_mass_fractions {csat_peridotite_poly_coeffs, csat_gabbro_poly_coeffs, \
147+
csat_MORB_poly_coeffs, csat_sediment_poly_coeffs
148+
};
149+
150+
std::vector<std::vector<double>> devolatilization_onset_temperatures {Td_peridotite_poly_coeffs, Td_gabbro_poly_coeffs, \
151+
Td_MORB_poly_coeffs, Td_sediment_poly_coeffs
152+
};
153+
};
154+
}
155+
156+
}
157+
}
158+
159+
#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)