# Perple_X 07 Seismic Velocity Calculations

!!!WARNING!!!

This tutorial is for an out-of-date version (Perple_X 07), the current version (Perple_X 666) of the tutorial is at seismic velocity tutorial.

## Contents

Introduction

Thermodynamic data file preparation

A comment on missing data (it's missing)

A Worked Example

Running Build: problem setup

Running Vertex: phase diagram section computation

Running Pssect: phase diagram section depiction

Running Werami: extraction of seismic velocities

Pscontor

Figures:

Figure 1: Perple_X program structure

Figure 2: Metabasalt phase relations

Figure 3: Metabasalt P-wave velocity

## Introduction

Seismic compressional and shear wave velocities through rock are determined by: the corresponding velocities within the constituent minerals; the amounts of the minerals; and rock texture. For a rock in thermodynamic equilibrium, the first two factors are determined entirely by thermodynamic properties and are computed automatically by Perple_X. The third factor, texture, determines the relation of the bulk rock seismic velocities to those of the constituent minerals. As texture is not a thermodynamic property, the averaging schemes used to compute bulk rock velocities are equivocal. In Perple_X, Voigt-Reuss-Hill (VRH) averaging is used whereby the aggregate velocity is computed from the average of the arithmetic and geometric means of the elastic moduli of the individual phases weighted by their volume fraction. The main advantage of the VRH scheme, which has no theoretical basis, is that it does not depend on texture. This averaging is done within the Perple_X program WERAMI.

Recognition of the thermodynamic constraints on seismic velocities has led some geophysicists to use free energy minimization to compute mineral modes and compositions for rocks. These modes and compositions are then used to estimate seismic velocities from semi-empirical models. This approach is perhaps justified by the uncertainties involved, but it is unquestionably thermodynamically inconsistent in that it amounts to using one thermodynamic database to establish the amounts of the minerals and another database to compute velocities. For geophysicists who remain nonetheless enamored with empirical estimates of elastic modulii, Perple_X can be used to make calculations in this manner. However, I personally advocate using the same thermodynamic data for the computation of both the stable mineralogy and seismic velocities. The approach actually used by Perple_X is determined by the structure of the data base as defined by the user and discussed below.

Typical thermodynamic data tabulations used for phase equilibrium calculations are designed for the computation of isostatic phase equilibria (i.e., for systems without differential stresses), however seismic wave propagation involves both shear and isostatic stresses. Therefore to compute seismic velocities with Perple_X conventional thermodynamic data bases must be expanded to define the shear modulus (m), which characterizes the thermodynamic response of each phase to shear stress. This document describes the necessary modifications to the thermodynamic data file and illustrates an application in which seismic velocities are recovered from a phase equilibrium calculation.

## Modification of the Thermodynamic Data File

This section is at: http://www.perplex.ethz.ch/perplex_seismic_velocity_I.html

## Computational Algorithm for Seismic Velocities

This section is at: http://www.perplex.ethz.ch/perplex_seismic_velocity_I.html

## A Worked Example

The example provided here (originally from Connolly & Kerrick 2002) illustrates the computation of absolute and relative seismic velocities through an altered metabasalt at p-T conditions relevant to subduction zones at depths of 90-240 km. Relative velocities are, in this case, the seismic velocities in the metabasalt divided by the corresponding velocities in the ambient mantle, which is represented here by a pyrolite model composition. The type of calculation required for this example is commonly referred to as a phase diagram section or sometimes simply as a pseudosection.

The calculation and interpretation of a phase diagram section consists of five steps (Figure 1): the problem definition; the phase diagram section computation and depiction; computational assessment; digital interpretation; and graphical interpretation. Each step involves different programs and is described in the following sections:

Figure 1. Perple_X program structure for gridded minimization calculations (i.e., for phase diagram section and phase fractionation calculations).

## Running BUILD

The computation of relative velocities requires the calculation of phase relations for both the metabasalt of interest and the pyrolite mantle model. The BUILD dialog below details the problem definition for the metabasalt, the slightly simpler dialog for the pyrolite calculation is in the file dialog17b.txt.

The following script details the prompts from the program BUILD (normal font) and user responses (bold font). Explanatory comments are written in red. Additional explanation on some prompts is available via the (help) links.

c\jamie\Berple_X>build

NO is the default ([cr]) answer to all Y/N prompts (help)

Enter name of problem definition file to be created, <100 characters, left justified [default = in]: (help)

in17a.dat

Enter thermodynamic data file name, left justified, [default = hp02ver.dat]: (help)

The file hp02ver.dat is a version of the Holland and Powell (1998, revised 2002) data base modified (Connolly 2005) to include shear modulii relevant to the high pressure mineralogy of metabasalts. The hp02ver.dat file is suitable for calculations to pressures approaching 12 GPa. The sfo05ver.dat data base, which is primarily from Stixrude & Lithgow-Bertelloni (2005, Khan et al. 2006), can be used over the entire range of mantle conditions, but is for a simplistic chemical model (see example #23).

`hp02ver.dat`
```
Enter the computational option file name, left justified, [default = perplex_option.dat]:
See: www.perplex.ethz.ch/perplex_options.html
```

Computational options used by Perple_X can be controlled via the computational option file. This file is normally named perplex_option.dat. Provision is made here to specify an alternative name because it sometimes desirable to have customized option files for specific purposes (e.g., high resolution or solvus calculations).

` `
`perplex_option.dat`
` `

The first major block of prompts from BUILD determine the manner in which the chemical components defined within the thermodynamic data file are to be used. VERTEX permits components to be treated in four different ways: saturated phase components are components whose chemical potentials are determined by the (assumed) stability of a phase, most commonly a fluid, containing these components; saturated components are components whose chemical potentials are determined by the (assumed) stability of a pure phase(s) and or solutions containing only the saturated-phase and saturated components; mobile components are components whose chemical potentials are defined explicitly; and thermodynamic components are components whose chemical potentials are the dependent variables of a phase diagram calculation. Phase diagram calculations require the specification of at least one thermodynamic component.

The current data base components are:

NA2O  MGO   AL2O3 SIO2  K2O   CAO   TIO2  MNO   FEO   O2    H2O   CO2

Transform them (Y/N)? (help)

n

Calculations with a saturated phase (Y/N)? (help)

The phase is: FLUID

Its compositional variable is: Y(CO2), X(O), etc.

` `

Here the system is not assumed to be saturated with respect to a fluid phase. With this closed system model fluid stability and composition is a function of the composition and physical conditions. This model requires that the volatile components of the fluid be specified as thermodynamic components as done below. There are problems with such a model, but the closed system model for water-CO2 bearing systems is certainly less arbitrary (and probably more realistic) than the petrological artifice of assuming that the system is saturated with respect to a fluid with a fixed composition. More realistic models for devolatilization require path-dependent models, such models can be implemented in Perple_X in “phase fractionation calculation” mode (and for phase fractionation with mass transfer see Connolly 2005). The main problem with the closed system model is that it implies that fluid generated by devolatilization remains within the system. This problem can be overcome by instructing WERAMI to compute the seismic properties of the solid aggregate (a limiting case).

n

Calculations with saturated components (Y/N)? (help)

n

Use chemical potentials, activities or fugacities as independent variables (Y/N)? (help)

n

Select thermodynamic components from the set: (help)

NA2O  MGO   AL2O3 SIO2  K2O   CAO   TIO2  MNO   FEO   O2    H2O   CO2

Enter names, left justified, 1 per line, [cr] to finish:

SIO2

TIO2

AL2O3

FEO

MGO

CAO

NA2O

K2O

H2O

CO2

**warning ver016** you are going to treat a saturated (fluid) phase component as a thermodynamic component, this may not be what you want to do. (help)

Select fluid equation of state: (help)

0 - X(CO2) Modified Redlich-Kwong (MRK/DeSantis/Holloway)

1 - X(CO2) Kerrick & Jacobs 1981 (HSMRK)

2 - X(CO2) Hybrid MRK/HSMRK

3 - X(CO2) Saxena & Fei 1987 pseudo-virial expansion

4 - Bottinga &Richet 1981 (CO2 RK)

5 - X(CO2) Holland &Powell 1991, 1998 (CORK)

6 - X(CO2) Hybrid Haar et al 1979/CORK (TRKMRK)

7 - f(O2/CO2)-f(S2) Graphite buffered COHS MRK fluid

8 - f(O2/CO2)-f(S2) Graphite buffered COHS hybrid-EoS fluid

9 - Max X(H2O) GCOH fluid Cesare & Connolly 1993

10 - X(O) GCOH-fluid hybrid-EoS Connolly & Cesare 1993

11 - X(O) GCOH-fluid MRK Connolly & Cesare 1993

12 - X(O)-f(S2) GCOHS-fluid hybrid-EoS Connolly & Cesare 1993

13 - X(H2) H2-H2O hybrid-EoS

14 - EoS Birch & Feeblebop (1993)

15 - X(H2) low T H2-H2O hybrid-EoS

16 - X(O) H-O HSMRK/MRK hybrid-EoS

17 - X(O) H-O-S HSMRK/MRK hybrid-EoS

18 - X(CO2) Delany/HSMRK/MRK hybrid-EoS, for P >10 kb

19 - X(O)-X(S) COHS hybrid-EoS Connolly & Cesare 1993

20 - X(O)-X(C) COHS hybrid-EoS Connolly & Cesare 1993

21 - X(CO2) Halbach &Chatterjee 1982, P >10 kb, hybrid-Eos

22 - X(CO2) DHCORK, hybrid-Eos

23 - Toop-Samis Silicate Melt

24 - f(O2/CO2)-N/C Graphite saturated COHN MRK fluid
25 - H2O-CO2-NaCl Aranovich & Haefner 2004

5

The second major block of prompts defines the explicit computational variables.

Specify computational mode: (help)

1 - Unconstrained minimization
2 - Constrained minimization on a 2d grid [default]
3 - Constrained minimization on a 1d grid
4 - Output pseudocompound data
5 - Phase fractionation calculations

Unconstrained optimization should be used for the calculation of composition, mixed variable, and Schreinemakers diagrams. Gridded minimization can be used to construct one- and two-dimensional phase diagram sections.

2

Here we are concerned with a normal P-T section; however, in gridded minimization it is possible to construct sections as a function of any arbitrarily defined composition (variable X(C1) in this prompt) or along a geothermal gradient, for examples of such calculations refer to the T-X, P-X and X-X section and phase abundances along an abritary path tutorials.

The data base has P(bars) and T(K) as default independent potentials. Make one dependent on the other, e.g., as along a geothermal gradient (y/n)? (help)

n

Select x-axis variable: (help)

1 - P(bars)

2 - T(K)

3 - Composition X(C1)* (user defined)

*X(C1) can not be selected as the y-axis variable

2

Enter minimum and maximum values, respectively, for: T(K)

673 1473

Enter minimum and maximum values, respectively, for: P(bars)

30000 80000

The following blurb provides information on the resolution with which phase relations are mapped as a function of the section variables. This resolution is controlled via the grid_parameters keywords specified in the perplex_option.dat file. When the auto_refine option is set on (the default) calculations are done in two iterations: an initial low quality "exploratory stage" is used to establish the range of phase compositions; and a second high resolution "auto-refine stage" to refine the final results. The two stages are executed in a single calculation if the auto_refine keyword is set to "auto".

For gridded minimization, grid resolution is determined by the number of levels (grid_levels) and the resolution at the lowest level in the X- and Y-directions (x_nodes and y_nodes) these parameters are currently set for the exploratory and autorefine cycles as follows:

stage        grid_levels  xnodes  ynodes    effective resolution
exploratory       1          20      20        20 x  20 nodes
auto-refine       4          60      60       473 x 473 nodes

To change these options edit or create the file perplex_option.dat
See: www.perplex.ethz.ch/perplex_options.html

Specify component amounts by weight (Y/N)?

y

Enter weight amounts of the components: (help)

SIO2  TIO2  AL2O3 FEO   MGO   CAO   NA2O  K2O   H2O   CO2

for the bulk composition of interest:

`45 1.1 15.26 9.84 6.54 12.66 2.03 0.55 2.63 2.9`
` `

The next block of prompts define the output required by the user.

Do you want a print file (Y/N)?

n

Enter the plot file name, <100 characters, left justified [default = pl]: (help)

plot17a

The final block of prompts from BUILD define the phases for the calculation.

Exclude phases (Y/N)? (help)

y

Here, experience leads me to make four exclusions: Clinozoisite (cz) is excluded to suppress the zoisite-clinozoisite transition, which is primarily a distraction in terms of interpreting phase equilibria. The water endmember of the Holland & Powell melt model (h2oL), is excluded because the endmember is (incorrectly) more stable than water at high pressure and low temperature. The "aluminum-free chlorite" endmember (afchl) of the Baker, Holland & Powell chlorite model is excluded because the use of this endmember requires many more pseudocompounds to resolve chlorite compositions yet it has almost no effect on the computed phase relations. The ordered-dolomite (odo) endmember required for Holland & Powell's (2004) dolomite model is excluded because the model, of dubious value in any case, is not employed here.

Do you want to be prompted for phases (Y/N)?

n

Enter names, left justified, 1 per line, [cr] to finish:

cz

afchl

h2oL

odo

Do you want to treat solution phases (Y/N)?

y

Enter solution model file name [default = solut_07.dat] left justified, <100 characters: (help)

`solut_07.dat`
` `
`...diagnostics omitted... `
` `

BUILD reads the solution model file and determines which solutions are consistent with the specified problem. The diagnostics (omitted here) can usually be ignored. They are useful if you discover that a solution you wish to include is not present in the list of solutions BUILD offers.

Select phases from the following list, enter 1 per line, left justified, [cr] to finish: (help)

Bio(HP)     TiBio(WPH)  TiBio(HP)   Chl(HP)     O(SG)       E(SG)

Omph(HP)    melt(HP)    pMELTS(G)   Pheng(HP)   Sapp(HP)    Sapp(KWP)

Osm(HP)     GaHcSp      F           F(salt)     T           St(HP)

Ctd(HP)     Carp        hCrd        Sud(Livi)   Sud         Cumm

Anth        Gl          Tr          TrTsPg(HP)  GlTrTs      GlTrTsPg

Amph(DHP)   Amph(DPW)   feldspar    Pl(h)       Kf          San

MaPa        KN-Phen     MuPa        Mica(CH)    Bio(TCC)    O(HP)

Cpx(l)      Cpx(h)      Mont        Do(HP)      M(HP)       Do(AE)

Cc(AE)      Sp(JR)      Sp(GS)      Sp(HP)      IlGkPy      Neph(FB)

GrPyAlSp(B  GrPyAlSp(G  Gt(HP)      Gt(EWHP)    Gt(WPH)     A-phase

Chum        Atg         B           P           Qpx         Mn-Opx

OrFsp(C1)   AbFsp(C1)   Pl(I1,HP)   Fsp(C1)     oCcM(HP)    O(stx)

Sp(stx)     Gt(stx)     Opx(stx)    Cpx(stx)    o-Amph      cOmph(HP)

CrGt        CrOpx(HP)   Opx(HP)     Cpx(HP)     CrSp        Ca-Amph(D)

Na-Amph(D)

F

Chl(HP)

Mica(CH)

San

Pl(h)

Gt(HP)

Do(HP)

M(HP)

Opx(HP)

Amph(DPW)

Omph(HP)

Bio(TCC)

T

Because fluid is not specified as a saturated phase it is necessary to specify explicitly that a fluid solution may be present (i.e., the solution model F).

Enter calculation title:

staudigel closed system model

Top

## Running VERTEX

If the auto_refine option is on calculations are done in two stages: an initial low quality exploratory stage establishes the range of phase compositions; and a final high resolution auto-refine stage. These stages are executed in a single calculation, as done here, only if the value of the auto_refine keyword is auto.

Before running VERTEX be sure to verify that the auto_refine keyword is set to auto in your copy of perplex_option.dat.

c\jamie\Berple_X> vertex

VERTEX begins by echoing the current settings for computational options, these are either the defaults or the values specified in perplex_option.dat.

Perple_X options are currently set as:

Keyword:               Value:     Permitted values [default]:
dependent_potentials   off        off [on ]
solvus_tolerance       0.050      0->1 [0.05], 0 => p=c pseudocompounds, 1 => homogenize
zero_mode              0.1E-05    0->1 [1e-6], < 0 => off
zero_bulk              0.1E-05    0->1 [1e-6], < 0 => off
iteration limit      4         >0 [3], value 1 of iteration keyword
refinement factor    3         >2*r [5], value 2 of iteration keyword
refinement points   10         0->100 [10], value 3 of iteration keyword
initial_resolution     0.10       0->1 [0.1], 0 => off
stretch_factor         0.016      >0 [0.0164]
reach_factor           0.9        >0.5 [0.9]
subdivision_override   off        [off] lin str
auto_refine            aut        off [manual] auto
auto_refine_factor_I    3         >=1 [3]
auto_refine_factor_II   5         >=1 [10]
auto_refine_factor_III  2         >=1 [3]
auto_refine_slop       0.010      [0.01]
hard_limits            off        [off] on
T_stop (K)                0.0     [0]

Gridded minimization parameters:
x_nodes              20 / 60    [20/40], >0, <2048; effective x-resolution   20 / 473 nodes
y_nodes              20 / 60    [20/40], >0, <2048; effective y-resolution   20 / 473 nodes
grid_levels           1 / 4     [1/4], >0, <10
1d_path              10 /150    [20/150], >0, <2048
Parameters for Schreinemakers and Mixed-variable diagram calculations:
variance              1 /99     [1/99], >0, maximum true variance
increment         0.100 /0.025  [0.1/0.025], default search/trace variable increment
efficiency             3        [3] >0 < 6
reaction_format      min        [min] full stoichiometry S+V everything
reaction_list        off        [off] on
console_messages     on         [on] off
short_print_file     on         [on] off
WERAMI/MEEMUM output options:
compositions           mol        wt  [mol]
proportions            vol        wt  [vol] mol
interpolation          on         off [on ]
extrapolation          off        on  [off]
vrh_weighting          0.5        0->1 [0.5]
cumulative_modes       off        [off] on

Worst case (Cartesian) compositional resolution (mol %)
Exploratory     0.146       10.000          10.000
Auto-refine     0.049        2.000           5.000
* - composition/mixed variable diagrams, ** - Schreinemakers diagrams

To change these options edit or create the file perplex_option.dat
See: www.perplex.ethz.ch/perplex_options.html

VERTEX then requests the problem definition file created by BUILD (in17a.dat):

Enter problem definition file name (i.e. the file created with BUILD), left justified:

in17a.dat

and lists the input/output files to be used/created for the calculation. Several output files are created in addition to the plot file "plot17a" requested by the user, among these are: "bplot17a" named by adding the prefix "b" to the plot file name, "pseudocompound_glossary.dat", and "auto_refine_in17a.dat.txt". The files pseudocompound_glossary.dat and auto_refine_in17a.dat.txt (named by appending the name of the problem definition file to the prefix "auto_refine_" and adding the suffix ".txt") contain information about the range of compositions considered in the calculation. In addition to the auto_refine text file, VERTEX generates a less user friendly version of the same data in "auto_refine_in17a.dat" (named by appending the problem definition file name to the prefix "auto_refine_"), which is used by VERTEX in the auto-refine stage. It is only necessary to be aware of these files if you wish to store or analyze results in a different directory than the working directory for VERTEX.

Reading thermodynamic data from file: hp02ver.dat

Writing print output to file: none requested

Writing plot output to file: plot17a

Writing pseudocompound glossary to file: pseudocompound_glossary.dat

Reading solution models from file: solut_07.dat

Writing bulk composition plot output to file: bplot17a

Writing data for auto-refinement to file: auto_refine_in17a.dat.txt

` `

All "make definitions" consistent with the problem definition file are output next. Make definitions define thermodynamic entities (e.g., buffers, endmembers or reactions) as a linear combination of independent entries in the thermodynamic data file. The make definitions are in the header section of the thermodynamic data file (hp02ver.dat) and can be modified, deleted or created for your own purposes.

Summary of valid make definitions:

MGO   AL2O3 SIO2  K2O   TIO2  FEO   H2O
tbi       1.00  0.50  1.00  0.50  1.00  0.00  1.00
tbi1      2.00  0.50  3.00  0.50  1.00  0.00  0.00
fbr       0.00  0.00  0.00  0.00  0.00  1.00  1.00
fchum     0.00  0.00  4.00  0.00  0.00  9.00  1.00
fphA      0.00  0.00  2.00  0.00  0.00  7.00  3.00
fper      0.00  0.00  0.00  0.00  0.00  1.00  0.00
fatg      0.00  0.00 34.00  0.00  0.00 48.00 31.00
sil8L     0.00  1.60  1.60  0.00  0.00  0.00  0.00
fo8L      4.00  0.00  2.00  0.00  0.00  0.00  0.00
fa8L      0.00  0.00  2.00  0.00  0.00  4.00  0.00
q8L       0.00  0.00  4.00  0.00  0.00  0.00  0.00

The remaining section of the initial console output (abridged below) consists of information and warnings on the solution models to be used in the calculation. These messages inform the user as to which solution models in the solution model file are compatible with the problem definition, i.e., that the solution end-members are possible phases for the calculation. If only a subset of the solution model end-members are present, the messages indicate how the model will be reformulated in terms of this subset. If you do not get the solution models you requested, or if you get different solution compositions than you expect, these messages may be helpful. The diagnostics also indicate the number of static pseudocompounds used for each solution, this number is useful for determining which solution models are costly in terms of computer memory.

**warning ver114** the following endmembers are missing for Bio(TCC)
ffbi_i    fbi

**warning ver050** reformulating reciprocal solution: Bio(TCC) because of missing endmembers.(reformulation can be controlled explicitly by excluding additional endmembers).

338 pseudocompounds generated for: Bio(TCC)

**warning ver114** the following endmembers are missing for Chl(HP)
mame_i    mafchl_i  mnchl     fafchl_i  afchl

**warning ver050** reformulating reciprocal solution: Chl(HP)    because of missing endmembers. (reformulation can be controlled explicitly by excluding additional endmembers).

118 pseudocompounds generated for: Chl(HP)
63 pseudocompounds generated for: Omph(HP)
9 pseudocompounds generated for: F
118 pseudocompounds generated for: T

**warning ver114** the following endmembers are missing for Amph(DPW)
mfets     ffets_i

**warning ver050** reformulating reciprocal solution: Amph(DPW)  because of missing endmembers. (reformulation can be controlled explicitly by excluding additional endmembers).

3141 pseudocompounds generated for: Amph(DPW)
9 pseudocompounds generated for: Pl(h)
9 pseudocompounds generated for: San
4288 pseudocompounds generated for: Mica(CH)
9 pseudocompounds generated for: Do(HP)

**warning ver114** the following endmembers are missing for M(HP)
rhc
9 pseudocompounds generated for: M(HP)

**warning ver114** the following endmembers are missing for Gt(HP)
spss
63 pseudocompounds generated for: Gt(HP)
118 pseudocompounds generated for: Opx(HP)

Total number of pseudocompounds:    8292

** Starting exploratory computational stage **

The preceding line terminates the initial section of the console output. Often this section is followed by warning messages such as:

**warning ver369** failed to converge at T= 673.00 K, P= 72105.3 bar Using Murnaghan EoS, probably for Ghiorso et al. MELTS/PMELTS endmember data. Volume estimated using 3rd order Taylor series.

that indicate unstable numerical behavior of specialized equations of state, most commonly these result because the calibration is being used beyond the range of physical conditions for which it was intended. In general, these warnings have no significant consequences, but if you have doubts it is possible to eliminate the offending entity (e.g., here endmembers of the pMELTS(G) silicate melt model are causing problems at subsolidus temperatures and pressures far beyond the 4 GPa upper limit for the pMELTS calibration).

During the computation VERTEX outputs progress information:

5.0% done with low level grid.
10.0% done with low level grid.
15.0% done with low level grid.
...
85.0% done with low level grid.
90.0% done with low level grid.
95.0% done with low level grid.
100.0% done with low level grid.

Beginning grid refinement stage.

After completion of the exploratory stage, VERTEX indicates which solution phases were not stable and the compositional ranges of the stable solution phases, this information is echoed in the "auto_refine_in17a.dat.txt" file and is used to increase for the subsequent auto-refine stage.

WARNING: The following solutions were input, but are not stable:

Bio(TCC)
Chl(HP)
Amph(DPW)
Pl(h)
Opx(HP)

Endmember compositional ranges for model: Omph(HP)

Endmember   Minimum   Maximum
di          0.40153   0.82889
jd          0.07778   0.52822

Endmember compositional ranges for model: F

Endmember   Minimum   Maximum
CO2         0.00000   0.33032

Site fraction ranges for multisite model: T

Site 1
Species   Minimum   Maximum
1      0.91556   0.94222

Site 2
Species   Minimum   Maximum
1      0.99444   0.99889

...output abridged...

At this point if the auto_refine keyword were set to manual (i.e., "man"), VERTEX would terminate, allowing the user to examine intermediate results, but requiring that VERTEX be run a second time to obtain the final high resolution calculation. However this example has been run with auto_refine set to auto (i.e., "aut"), thus the VERTEX automatically begins the auto-refine stage of the calculation:

Reading auto-refinement data from file: auto_refine_in17a.dat

At the beginning of auto-refine stage VERTEX eliminates all solution phases specified in the problem definition file (in17a.dat) that were not stable in the exploratory calculation and then restricts the compositional ranges for the remaining solutions to those found in the exploratory calculation (plus an increment [auto_refine_slop] to allow for inaccuracy), additionally the compositional resolution of the calculation in these restricted compositional resolution is increased by the auto_refine_factor.

Eliminating solution model: Chl(HP)    in auto-refinement.

Eliminating solution model: Pl(h)      in auto-refinement.

Eliminating solution model: Opx(HP)    in auto-refinement.

Eliminating solution model: Amph(DPW)  in auto-refinement.

Eliminating solution model: Bio(TCC)   in auto-refinement.

148 pseudocompounds generated for: Omph(HP)
11 pseudocompounds generated for: F
2 pseudocompounds generated for: T
31 pseudocompounds generated for: San
164205 pseudocompounds generated for: Mica(CH)
4 pseudocompounds generated for: Do(HP)
14 pseudocompounds generated for: M(HP)
117 pseudocompounds generated for: Gt(HP)

Total number of pseudocompounds:  164532

** Starting auto_refine computational stage **

During auto-refinement VERTEX may write diagnostics such as:

WARNING: composition of solution Do(HP) has reached an internal limit (0.759) on site 1 for species 1.

If this warning occurs during auto-refinement, the problem can be circumvented by increasing the auto_refine_slop specified in perplex_option.dat.

For the remainder of this calculation the internal limit has been relaxed to: 0.726

that commonly indicates that the composition of a phase has exceeded the range discovered in a previous exploratory stage. When this occurs VERTEX relaxes the range as indicated, which usually prevents any significant consequences. If this diagnostic occurs during the exploratory stage itself it may have more serious consequences. In this case, it is wise to widen the relevant limits specified in the solution model file.

...

98.3% done with low level grid.

100.0% done with low level grid.

In contrast to the exploratory stage, the grid_parameter keyword specifies a multi-level (i.e., 4-level) grid for the auto-refine stage, thus after the calculation of the first (or "low level") grid, VERTEX makes three refinements of the grid, each of these refinements takes roughly twice the time of the previous refinement.

Beginning grid refinement stage.

649 grid cells to be refined at grid level 2

...working (   501 minimizations done)

...working (  1003 minimizations done)

refinement at level 2 involved   1376 minimizations

4976 minimizations required of the theoretical limit of  14161

1202 grid cells to be refined at grid level 3

...working (   128 minimizations done)

...working (   629 minimizations done)

...working (  1131 minimizations done)

...working (  1634 minimizations done)
...

VERTEX terminates as in the exploratory cycle with a summary of the phase compositions and writes a new auto-refine data file.

## Running PSSECT

PSSECT generates PostScript graphical output and can be used once the phase diagram section has been calculated with VERTEX. The default plot from PSSECT shows the computed phase relations (Figure 3) and is generated by the following dialog.

c\jamie\Berple_X> pssect

...

Enter problem definition file name (i.e. the file created with BUILD), left justified:

in17a.dat

PostScript will be written to file:

plot17a.ps

Answer yes to the following prompt to modify features such as axis labeling, size, and shape.

Modify the default plot (y/n)?

n

` `
` `
` `
` `
` `
` `
` `
` `
` `
` `
` `
` `

Figure 2 Computed metabasalt phase relations (plot17a.ps), this figure differs from Figure 2a of Connolly & Kerrick (2002) because the published figure was calculated with different solution models.

## Running WERAMI

WERAMI permits the user to extract information from a phase diagram section at a point, throughout a region, or along a line or curve. The latter two options generate output that can be plotted with programs PSCONTOR or PSPTS or imported into general-purpose programs such as Matlab. The dialog and output for each of these modes is detailed in the . Here the second mode is used to extract P-wave velocities for the metabasalt over the entire range of pressure-temperature conditions represented by the phase diagram section.

c\jamie\Berple_X> werami

Enter problem definition file name (i.e. the file created with BUILD), left justified:

in17a.dat

Select operational mode:

1 - compute properties at specified conditions

2 - create a property grid (plot with pscontor)

3 - compute properties along a curve or line (plot with pspts)

4 - as in 3, but input from file

0 - EXIT

2

Select a property:

1 - Specific enthalpy (J/m3)

2 - Density (kg/m3)

3 - Specific Heat capacity (J/K/m3)

4 - Expansivity (1/K, for volume)

5 - Compressibility (1/bar, for volume)

6 - Weight percent of a component

7 - Mode (Vol %) of a compound or solution

8 - Composition of a solution

9 - Grueneisen thermal ratio

11 - Shear modulus (bar)

12 - Sound velocity (km/s)

13 - P-wave velocity (Vp, km/s)

14 - S-wave velocity (Vs, km/s)

15 - Vp/Vs

16 - Specific Entropy (J/K/m3)

17 - Entropy (J/K/kg)

18 - Enthalpy (J/kg)

19 - Heat Capacity (J/K/kg)

20 - Specific mass of a phase (kg/m3-solid)

21 - Poisson ratio

22 - Molar Volume (J/bar)

23 - Chemical potentials (J/mol)
24 - Assemblage Index
25 - Modes of all phases

13

Calculate individual phase properties (y/n)?

n

Include fluid in computation of aggregate (or modal) properties (y/n)?

n

Change default variable range (y/n)?

n

Enter number of nodes in and x and y directions: (help)

100 100

Writing grid data to file: cplot17a

**warning ver637** Immiscibility occurs in one or more phases interpolation will be turned off at all affected nodes. To overide this feature at the risk of computing inconsistent properties turn solvus testing off in perplex_option.dat

Range of >0 data is: 7.90399 -> 8.88201

n

Select operational mode:
1 - compute properties at specified conditions
2 - create a property grid (plot with pscontor)
3 - compute properties along a curve or line (plot with pspts)
4 - as in 3, but input from file
0 - EXIT

0

The structure of the grid data file is described in README WERAMI GRID FORMAT.

## Graphical Representation

Grid data generated by WERAMI can be plotted by various commercial packages such as MatLab or Maple, alternatively (the relatively crude) Perple_X program PSCONTOR can be used to generate an interpreted PostScript plot. Both possibilities are demonstrated here.

### PSCONTOR

c\jamie\Berple_X>pscontor

Enter the CONTOUR plot file name:

cplot17a

PostScript will be written to file:

cplot17a.ps

Reset plot limits (y/n)?

n

Modify default drafting options (y/n)?

- picture transformation

- x-y plotting limits

- relative lengths of axes

- text label font size

- line thickness

- curve smoothing

- axes labeling and gridding

y

Modify default picture transformation (y/n)?

n

Modify ratio of x to y axis length (y/n)

n

Modify default text scaling (y/n)?

n

Use minimum width lines (y/n)?

n

Fit curves with B-splines (y/n)?

y

Contoured variable range: 7.940270 ->8.953970

Modify default contour interval (y/n)?

y

Enter min, max and interval for contours:

7.95 8.95 0.05

Modify default axes (y/n)?

y

Enter the starting value and interval for major tick marks on

the X-axis (current values are: 673.0 160.0)

Enter new values:

673 160

Enter the starting value and interval for major tick marks on

the Y-axis (current values are: 0.300E+05 0.100E+05)

Enter new values:

30000 10000

Cancel half interval tick marks (y/n)?

n

Draw tenth interval tick marks (y/n)?

n

Rescale axes text (y/n)?

n

Turn gridding on (y/n)?

n

The output (Figure 4) from PSCONTOR can be printed on any PostScript printer or viewed on a Postscript viewer such as Ghostview. Alternatively, the postscript plot may be imported into graphical editors such as Adobe Illustrator or Coreldraw. PSCONTOR does not label the plotted contours, but it does use a heavy solid line for the minimum contour and a heavy broken line to indicate the maximum contour value. If this information is not sufficient to determine the contour values, WERAMI can be used in mode 1 to query the value of the property of interest at any point within the phase diagram section.

### MatLab Plotting

The grid data file generated by WERAMI (README WERAMI GRID FORMAT) can be imported into MatLab or other commercial plotting routines. Three MatLab plotting scripts (matlab_plotting_scripts.zip) are provided with Perple_X to allow users to view grid data quickly:

surface_plot.m                         - plot a 3D image of the data

automatically_labeled_contour_plot.m   - 2D contour plot, contours chosen and labeled by Matlab

manually_labeled_contour_plot.m        - 2D contour plot, contours chosen and labeled by user

ratio_contour_plot_manually_labeled.m  - 2D contour plot of the ratio of properties in two different

WERAMI plot files (e.g., density ratio of two rocks)

Examples of the first and third plots are shown in Figure 4. Users unfamiliar with MatLab should be aware that MatLab expects the user to interactively label the contours generated by the third script by "clicking" on the contours, once the user has finished she must terminate the interactive labeling mode by pressing enter.

To make a MatLab plot, the user must start MatLab and change directory to the directory containing the Perple_X plotting scripts and then type the name of the desired script (without the “.m” suffix), e.g., “surface_plot”, the rest is self-explanatory.

Figure 3: Three plot variations of the computed metabasalt P-wave velocities (km/s). The uppermost plot is an interactively labeled contour plot generated in MATLAB by the script manually_labeled_contour_plot. The middle plot is a false colour image also generated with surface_plot.m, but with the lighting changed from “phong” to “giraud”. The lowermost plot is the plot generated by PSCONTOR.