Perple_X Equilibrium Reactive Transport Calculations: FRAC2D Model Configuration
Thoughts on configuring/running a FRAC2D calculation
Explanation/Definition of Variables
Computation/Plotting of Mass Fluxes
Figure 1: FRAC2D vs TITRATE reactive transport models
Figure 2: Elemental fluxes across the slab-mantle interface
Figure 4: Total solute-, CH4-, and CO2-content of the fluid
Figure 5: Bulk column sulfur content
Perple_X is configured to do two general types of equilibrium reactive transport models, FRAC2D and TITRATE (Figure 1). Both configurations are referred to within Perple_X as 2d-phase fractionation models (Connolly 2005). In these calculations a 1-dimensional column of rock is subject to perturbations that lead to the formation of one or more mobile phases (e.g., melt and low-density electrolytic or molecular fluids). The rate at which the perturbations are applied adds a second effective dimension, time, to the model. The mobile phase(s), i.e., fractionated phase, moves upward through the column until it is consumed or passes through the top of the column. The two types of calculations are distinguished by the mechanism for perturbing the initial equilibrium condition: in FRAC2D calculations, the perturbation is caused by changing pressure-temperature conditions, as occurs within a subducted slab or during regional metamorphism; in TITRATE calculations the perturbation is a mass flux (e.g., the presumed composition of a mobile phase) to the base of the column, as occurs in the mantle wedge during subduction and other metasomatic processes. This page documents the FRAC2D configuration, refer to the TITRATE page for an outline of TITRATE configuration. FRAC2D and TITRATE represent limiting cases, the intermediate case in which the column moves through a dynamic P-T field and is simultaneously subject to a basal mass flux can be computed as a variant of the FRAC2D configuration in which the TITRATE variable (see below for definition of all input/output variables) is set.
The input files for the FRAC2D example outlined here are at www.perplex.ethz.ch/perplex/examples/FRAC2D are designated by the project name frac2d. The fractionated, or mobile phase, is represented by the solution model COH-Fluid+, a model for low-density C-O-H-S fluids. Depending on the aq_lagged_speciation option specified in the Perple_X option file, COH-Fluid+ may represent either molecular or electrolytic fluids (Connolly and Galvez, 2018). The example here uses the molecular version of the model and thus corresponds to an up-date of the subduction model used by Connolly (2005). It is advisable to work with the molecular model before attempting the electrolytic model which introduces additional complexity and uncertainty. As in all Perple_X calculations, the project name has no special significance, the name frac2d used here is chosen solely for purpose of the example. In this example, equilibrium phase relations are computed for a column of rock representing the lithologies of the upper 24.5 km of an oceanic slab as it is subducted incrementally from 80 to 200 km depth, the column is oriented vertically because it is assumed this is the direction of fluid flow. Consequently, the height of the column is 24.5 km / cosine (45 degrees) ~ 34.6 km (see Fig 6 of Connolly 2005), where 45 degrees is the model slab dip. For purposes of the calculation, the column is discretized as a set of representative nodes. After each depth increment, beginning from the lowest representative node, equilibrium phase relations are computed. If the mobile phase is stable, its mass is removed (fractionated) from the node and added to the overlying node. The procedure is then repeated for the overlying node until, ultimately, the top of the column has been reached. This model carries strong assumptions about the nature of fluid generation and flow. Specifically, these assumptions are that fluid generation proceeds from the base of the column upwards, that fluid flow is instantaneous on the subduction time scale (DT0) and that the flow is pervasive and vertical. There is little basis for the first assumption, although the case can be made that any fluid accumulation at depth will induce hydro-mechanical instabilities that causes the fluid at depth to catch-up fluid at shallower depths (Connolly 1997). Experiments with alternative assumptions, e.g., that fluid generation proceeds from the top downward, show that this assumption does not have profound consequences for mass transport. Pervasive flow maximizes the efficacy of mass transport through the slab and is therefore a useful extremal scenario. The alternative extremal scenario, that fluid flow is perfectly channelized by fractures or tube-like porosity waves (Connolly & Podladchikov, 2007, 2013; Connolly 2010) corresponds to simple Rayleigh fractionation (e.g., Section 4.4 of Connolly & Galvez 2018) and is not considered further here. The independent variables for the calculation the depth of the slab/mantle interface (Z0) and the vertical height above the base of the slab column (DZ), the pressure-temperature field from a geodynamic model is expressed as a function of Z0 and DZ (see below, also Fig 7 of Connolly 2005). The model neglects the small effect of the divergence in the mass flux on density and, for the purposes of relating pressure and depth, assumes a constant density throughout the column and overlying rocks.
The primary input file for the calculation, frac2d.dat, is generated using BUILD, by choosing computational mode 7 (2-d Phase fractionation). In contrast to most Perple_X calculations in the case of FRAC2D calculations, this input file specifies only the data sources, thermodynamic components, and solution models to be considered during the calculation. A second auxiliary file, frac2d.aux, specifies the column lithologies, the model geometry, and the data necessary to define the pressure-temperature field during subduction. The auxiliary input file is generated manually, the primary purpose of this document is to explain the data in the frac2d.aux example so that it can be used as a template for new calculations.
Files at www.perplex.ethz.ch/perplex/examples/FRAC2D:
frac2d_noaq.dat - the primary input file (Figure 3).
frac2d_noaq.aux - the auxiliary input file (Figure 3).
frac2d_noaq_perplex_option.dat - the run time option file (Figure 3).
frac2d_aq.dat - variant with lagged-aqueous speciation enabled (Figures 2 and 4), the primary input file.
frac2d_aq.aux - variant with lagged-aqueous speciation enabled (Figures 2 and 4), the auxiliary input file.
frac2d_aq_perplex_option.dat - variant with lagged-aqueous speciation enabled (Figures 2 and 4), the run time option file.
solution_model.dat - the solution model file.
DEW13HP622ver_elements.dat - the thermodynamic data file.
NOTES: 1) The DEW13 aqueous solute data-base was not converted to Perple_X format using the default options specified in the DEW spreadsheet, it is used here only for purposes of demonstration, the DEW17 aqueous solute data-base should be used for quantitative applications. 2) These files can only be read by 6.8.6+ versions of Perple_X.
Steps in a FRAC2D calculation:
A) Run BUILD to create the primary input file, e.g., frac2d.dat.
B) Create the auxiliary input file, e.g., frac2d.aux.
C) Run VERTEX, answer yes to the fractionate phases prompt and fractionate the mobile phases of interest (e.g., COH-Fluid+).
D) Plot and/or analyze phase relations using PSSECT and WERAMI (see section 3, Figures 2, 3, 4). Note that in fractionation calculations, the data reported by VERTEX include the amount and composition of the mobile phase prior to fractionation. In FRAC2D calculations, the amount of the mobile present at a given height within the column is the product of the mobile phase generated at that height plus any mass gained by fractionation of the mobile phase at greater depths.
Thoughts on configuring/running a FRAC2D calculation:
E) If you are unfamiliar with Perple_X, familiarize yourself with the programs by following a tutorial such as www.perplex.ethz.ch/perplex_66_seismic_velocity.html.
F) Set auto_refine to off or manual, using the low-resolution (exploratory stage) result will identify problems quickly, these problems must be resolved before a high-resolution (auto-refine stage) result is attempted.
G) Before configuring your own calculation, begin at step C, above, using the input files at www.perplex.ethz.ch/perplex/examples/frac2d. Verify that you are able to obtain, understand, and analyze the output generated by the example.
H) Increase complexity gradually: Do not begin with an electrolytic multi-component fluid, rather begin with water, add molecular volatiles, and lastly add solutes. Do not begin with complex lithologic chemistry, rather begin with elements in a single oxidation state and exclude minor and/or poorly constrained elements. Construct phase diagram sections for each lithology within your model rock column over the entire pressure-temperature range of interest. If the phase relations in these sections do not make sense, then the phase relations in a FRAC2D calculation certainly will not make sense.
Explanation/Definition of FRAC2D Variables
Variables used to describe the FRAC2D calculation, unless otherwise indicated these are specified in the frac2d.aux file or computed on output:
AREA - the cross-sectional area of the subducted rock column, AREA = MASS/ZBOX/RHO.
ANNEAL - if .true., then at the initial condition of the calculation: a) if fluid is stable at any node in the column, it is removed and the composition of the node adjusted accordingly; b) the representative masses (MASS) of all nodes are renormalized to 1 kg.
C - the number of thermodynamic components specified in frac2d.dat
CLAY_1_1 ... CLAY_1_C
... ... ...
CLAY_N_1 ... CLAY_N_C - the molar compositions of the N column lithologies, from the base of the column upward. The components in the composition must be entered in the same order as the components are entered in frac2d.dat. For calculations in which the fluid is electrolytic it is advisable to use elemental components. See www.perplex.ethz.ch/perplex/faq/oxide_to_elemental_component conversion.pdf for an illustration of the oxide-to-elemental conversion used for the frac2d example. Note that these compositions are extensive. If ANNEAL is .true., the mass of the composition will be automatically adjusted to 1 kg. If ANNEAL is .false., the molar mass of each layer composition should be comparable.
C_1_0 ... C_1_NORD
... ... ...
C_NPOLY_0... C_NPOLY_NORD - the coefficients of the NPOLY NORD'th order polynomials that define the within slab geotherms, such that T(K) on the N'th geotherm is T_N = sum(C_N_I*ZI, I = 0,...,NORD). NOTE: Perple_X does (cannot) not test that the polynomials are valid, therefore it is advisable to plot the geotherms over the entire depth range relevant to a calculation. In some cases it may be impossible to represent a geotherm with a single polynomial, when this is the case use the PZFUNC option to specify the geotherms or explicit P(Z0,DZ)-T(Z0,DZ) functions.
DIP - slab dip in degrees
DPDZ - the pressure gradient with depth (bar/m). DPDZ = RHO(kg/m3)*g(m2/s)*(1 bar/105 Pascal), where g is gravitational acceleration (10 m/s2) and RHO is the average density of the column and any overlying rock. Perple_X can be modified to compute a thermodynamically consistent pressure for a loaded self-gravitating column; such complications seem unwarranted in light of other sources of uncertainty.
DT0 - the time interval (s) corresponding to one depth increment of the column, DT0 = DZ0/cosine(DIP)/V_SUB
DYNAMIC - if .true. Pressure and temperature are a function of the Z0-DZ coordinate. If DYNAMIC is .false., pressure and temperature are solely a function of DZ.
DZ - height above the base of the column; as DZ = 0 at Z0, DZ is negative within the column (see ZBOX).
DZ0 - the incremental change in depth of the column during subduction, DZ0 = (Z0MIN-Z0MAX)/NINT
DZGEO_1...DZGEO_NPOLY - the within slab (orthogonal) depths of the NPOLY geotherms that define the within slab P-T field.
MASS - the mass (g) associated with a representative node within the column. If ANNEAL is true, MASS is automatically adjusted to 1 kg; otherwise MASS is the mass corresponding to the extensive molar composition of the layer composition (CLAY_I_J, J = 1...C).
MASS_I - the mass (g) of a component transported by mobile phase I (e.g., fluid and/or melt).
NINT - the number of depth increments as the column is subducted from Z0MIN to Z0MAX. NINT = 1d_path_nodes - 1, where 1d_path_nodes is specified in frac2d_perplex_option.dat.
NORD - the order of the polynomials used to describe the within slab P-T field.
NPOLY - the number of geotherms used to describe the within slab P-T field, three should be adequate for most purposes. The within column temperature is computed by fitting an (NPOLY-1)'th order polynomial to the geotherms (see Connolly 2005).
RHO - the average density (kg/m3) of the column and overlying rock (see DPDZ), RHO = DPDZ/g*(105 Pa/bar)
P - pressure (bar) at any point within the column, P = Z*DPDZ.
Q - the flux (g/m2-s) of a component transported by mobile phase I through the column over time interval DT0 (e.g., fluid or melt) Q = MASS_I/AREA/DT0 (see step D below).
QINT - the time-integrated flux (g/m2) transported through the column by mobile phase I (see step D below).
PZFUNC - if .true., then internal functions are used to compute the P-T field of the column, e.g., from analytical functions (Davies, GJI, 1999) or tabulated results from, e.g., a geodynamic model. This option requires that the user program the functions.
T - temperature (K).
TITRATE - if .true., then the column is subject to a basal mass flux.
VERBOS - if .true., then VERTEX prints the amount and composition of the mobile phase at each node to the console. Additionally, if ANNEAL is also .true., then VERTEX outputs the initial composition at each node and the average composition of each layer after annealing.
V_SUB - the subduction rate (m/s).
V_HOR - the horizontal component of the slab velocity (m/s), V_HOR = V_SUB*cosine(DIP).
Z - the absolute depth (m) of any point within the column, as DZ increases upward, Z = Z0-DZ.
ZBOX - the vertical unit of discretization for the column in meters, i.e., the rock column is composed of boxes with a vertical dimension of ZBOX. The representative nodes of the column are located at DZ_I = -ZTOT+ZBOX*(1/2+[I-1]); thus the representative nodes at the top and bottom of the column are at DZ_TOP = -ZBOX/2 and DZ_BOT = -ZTOT+ZBOX/2, respectively. NOTE: the nodal coordinates of a layer should be entered in WERAMI or PSSECT to compute fluxes, e.g., to compute the flux across the top of the column, compute properties at column depth -ZBOX/2.
ZLAY_1...ZLAY_N - the vertical thicknesses (m) of the N lithological layers within the column, numbered from the base of the column upward. The vertical thickness is equal to the thickness of the layer orthogonal to the slab interface divided by the Cosine of the slab dip (DIP). NOTE that the layers are discretized with a resolution ZBOX, thus the layer thicknesses should be an integer multiple of ZBOX. If this is not true, then the model layer thickness will be (ZLAY_I\ZBOX)*ZBOX, where \ denotes the integer quotient. When a computation is run in VERTEX, the actual layer thicknesses and the nodal coordinates of the top and bottom of each layer are output to the console.
ZSURF - value of the depth coordinate at the surface (m), used for the computation of pressure at depth, i.e., P(bar) = (Z0-DZ)*DPDZ.
ZTOT - the total vertical thickness of column (m), ZTOT = sum(ZLAY_I, I = 1...N).
Z0 - depth (m) at the top of the column (m), i.e., the depth of the slab/wedge interface.
Z0MIN - the initial depth (m) of the slab/mantle interface
Z0MAX - the final depth (m) of the slab/mantle interface
Figure 2: Time-integrated elemental fluxes across the slab-mantle interface. The negligible fluxes of Fe, Al, and Mg are not shown. NOTES: 1) the fluxes must be scaled by ZBOX and RHO (see variable definitions above) to obtain rational units (e.g., g-mol/m2); 2) the option file linked from this page disables electrolytic fluid speciation so that carbon, sulfur, hydrogen, and oxygen are the only mobile elements. To reproduce the figure presented here the aq_lagged_speciation option must be set to T (true).
Computation and Plotting of Mass Fluxes
A) Set the absolute option to T in frac2d_perplex_option.dat.
B) Set the cumulative option to T in frac2d_perplex_option.dat if cumulative fluxes are desired, otherwise the fluxes are instantaneous, i.e., the average flux over the time interval DT0.
C) After setting the above options, start WERAMI:
i) Select operation mode 3 (properties a long a 1d path).
ii) Construct a linear profile, from (Z0 = Z0MIN, DZ = -ZBOX/2) to (Z0 = Z0MAX, DZ = -ZBOX/2).
iii) For optimal resolution specify NINT + 1 points for the profile.
iv) Select property 36 (all phase and/or system properties).
v) Specify option 2 (properties of a phase).
vi) Specify the name of the fractionated phase, e.g., COH-Fluid+.
vii) Answer y or n to the Include fluid in ... properties prompt, the answer is irrelevant for present purposes.
iix) The output is in the indicated tab file, e.g., frac2d_1.tab.
D) Plot the amount of the component of interest (e.g., S2_g) using perple_x_simple_plot (MATLAB), PSTABLE, or PYWERAMI. If the cumulative option is F or default, the plotted quantity corresponds to the MASS_I (g) of the component that is transported by the fluid across the slab/mantle interface (at DZ+ZBOX/2 = 0) as the top of the column is subducted from depth Z0 to Z0+DZ0. The instantaneous flux, i.e., the flux averaged over the time interval DT0, Q (g/m2-s) = MASS_I/AREA/DT0. If the cumulative option is T, then the plotted quantity is the cumulative mass TOT_MASS_I of the component transported by the mobile phase, the time-integrated mass flux QINT (g/m2) is TOT_MASS_I/AREA. The data in the tab file can be edited in a spreadsheet program so that Q and/or QINT are output directly. Alternatively, the axes of the plot may be rescaled to indicate the appropriate flux.
E) If a calculation is done with more than one mobile phase (e.g., melt and low-density fluid), step C must be repeated for each phase and the results summed to obtain the total mass flux transported by the mobile phases.
Figure 3: Phase relations within the subducted column as a function of column depth (Z0) and height within the column (DZ) plotted with PSSECT. Irregularities occur because small oscillations in the bulk composition at the representative nodes caused by metasomatism affect the stability of minor phases. These oscillations can be dampened by increasing NINT and decreasing ZBOX. The example shown here is for a calculation in which the fluid is solute-free (aq_lagged_speciation = F). Phase relations plotted directly with PSSECT for solute-bearing systems are often illegible, in such cases phase boundaries should be located by tabulating the modal abundances of the stable phase (WERAMI, mode 2, property choices 7, 25, or 36) and a plotting script (e.g., MATLAB or PYWERAMI).
<![if !vml]><![endif]><![if !vml]><![endif]><![if !vml]><![endif]>
Figure 4: Total solute molality and CH4- and CO2-content (mole fraction) in the fluid phase. Where no data is plotted the fluid is not stable. The data for this figure was generated with WERAMI (mode 2), property choice 40, and plotted in MATLAB with the perple_x_simple_plot script. The plot of methane concentration gives a false impression of the importance of such fluids in this example, i.e., the amount of the methane-rich fluid is minor.
Figure 5: Bulk sulfur content in the column. The pattern suggests that sulfur is scoured from the fluid as it infiltrates into progressively more silicic lithologies toward the slab surface. The three-dimensional perspective obscures the result that the maximum sulfur content is realized at the base of the sedimentary layer and that it subsequently decreases almost to the original sedimentary sulfur content at the slab interface. The data for this figure was generated with WERAMI (mode 2), property choice 6 (or 36) and plotted in MATLAB with the perple_x_simple_plot script.