Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A new approach to mass balance modeling: applications to igneous petrology Ke, Peiwan 1992

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_1993_spring_ke_peiwan.pdf [ 7.14MB ]
Metadata
JSON: 831-1.0052616.json
JSON-LD: 831-1.0052616-ld.json
RDF/XML (Pretty): 831-1.0052616-rdf.xml
RDF/JSON: 831-1.0052616-rdf.json
Turtle: 831-1.0052616-turtle.txt
N-Triples: 831-1.0052616-rdf-ntriples.txt
Original Record: 831-1.0052616-source.json
Full Text
831-1.0052616-fulltext.txt
Citation
831-1.0052616.ris

Full Text

A NEW APPROACH TO MASS BALANCE MODELING: APPLICATIONS TO IGNEOUS PETROLOGY  BY PEIWEN KE B.Sc., THE UNIVERSITY OF SCIENCE AND TECHNOLOGY OF CHINA, 1989  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTERS OF SCIENCE IN THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF GEOLOGICAL SCIENCES  We accept this thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA August 1992.  © Peiwen Ke, 1992.  In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission.  (Signature)  Department of  Geological Sciences  The University of British Columbia Vancouver, Canada  Date  DE-6 (2/88)  ^A v [ 2_6  , t  ,dy 3  ABSTRACT  An overdetermined linear system of equations is proposed to investigate mass balance relationships in igneous systems. A chi square (x 2 ) minimization technique is used to solve the set of linear equations. The x2 technique is superior to R 2 techniques because it considers the analytical errors on measurement data. Singular Value Decomposition (SVD) is used to search for the x 2 optimal solution. The attributes of the x 2 technique are discussed. The x2 minimization fitting technique, combined with SVD methods allows the confidence regions to be calculated and expressed by a geometric shape centered on the X2 optimal solution. The geometric shape varies with the dimension of the linear system: one-dimensional cases yield a line fragment; two-dimensional cases yield an ellipse and the three dimensions require an ellipsoid. The size of the confidence regions is determined by the confidence limits and is related to the coefficient matrix of the linear problem. Model with smaller confidence regions are considered to be better than models associated with larger confidence regions where the same confidence limit is used. Therefore, the confidence region can be used to select the model that best explains chemical differentiation in igneous system. Synthetic data sets have been analyzed using the new program MASSBAL.0 which solves the mass balance problems using the principles of the X 2 technique and SVD. The X2 technique presents a solution closer to the true one than does the traditional residual minimization fitting (R 2 technique). It overcomes the deficiencies that the R 2 technique has and appears to be a more robust and sensitive fitting technique.  II  The x 2 technique has been applied to the apparent clinopyroxene paradox formed in lavas from Diamond Craters, Oregon. Three models are proposed to explain the chemical difference: a three phase model (01+P1+Px), a two phase model (O1+Pl) and a one phase model (01). The X2 mass balance modeling results show that the three phase model (O1+Pl+Px) has the smallest minimum )(2 value and the most compact confidence region; thus it is preferred and clinopyroxene did participate in the fractionation process. The )( 2 technique has also been applied to study the assimilation problem in Paricutin volcano and the fractionation problem in PLA group of Mount Meager volcanic complex. The mass balance calculations show that assimilation has play a major role in Paricutin volcano and about 15% xenolith has contaminated the magma. In PLA fractionation processes, plagioclase and amphibole are shown to have caused the chemical difference between lavas, although other phenocrysts ( biotite and pyroxene ) did not participate.  III  TABLE OF CONTENTS  ABSTRACT ^  II  TABLE OF CONTENTS ^  IV  LIST OF TABLES ^  VI  LIST OF FIGURES ^  VIII  ACKNOWLEDGMENTS ^  X  CHAPTER 1 INTRODUCTION AND PRINCIPLES OF MASS BALANCE MODELING ^  1  1.1 Objectives ^  2  1.2 Previous Mass Balance Modeling Approaches ^  3  1.3 Fitting Techniques Used in Literature ^  9  CHAPTER 2 MASS BALANCE MODELING ALGORITHM ^ 11 2.1 Mass Balance Calculations in Igneous Differentiation ^ 11 2.2 Fractionation of Three Phases ^  12  2.3 Optimization of Mass Balance Equations ^  16  2.4 Mapped Confidence Region for the X 2 Optimization ^ 20 2.5 Goodness-of-Fit for Optimal Solutions ^  26  2.6 Computer Program and Calculation Procedures ^  26  CHAPTER 3 SIMULATED DATA SETS AND ALGORITHM TESTING ^ 31 3.1 Calculation of Simulated Data Set. ^  31  3.2 Simulated Data Sets with Associated Uncertainties. ^ 32 3.2.1 Solutions for Synthetic Data Sets with 1% Relative Error. ^ 34 3.2.2 The Effect of a on X 2 and R 2 Optimization ^ 39 3.2.3 x 2 Confidence Regions: A Measure of Goodness-of-Fit and A Test of Hypothesis ^ CHAPTER 4 NATURAL DATA SETS ^  IV  45 50  4.1 Paricutin Volcano: An Example of Assimilation ^ 51 4.2 Diamond Craters Volcano: An Example of Fractionation ^ 54 4.3 Mount Meager: Another Example of Fractionation ^ 57 4.2.1 Introduction ^  57  4.2.2 Mineralogy and Rock Chemistry Data of the Mount Meager Volcanic Complex ^  61  4.2.3 Mass Balance Calculations with the Plinth Assemblage (PLA) ^ 69 4.2.3.1 Comagmatic Relationships ^  69  4.2.3.2 Mass Balance Calculations ^  72  CHAPTER 5 SUMMARY AND DISCUSSION ^  78  REFERENCES ^  83  APPENDIX ^  88  Al Linear Algebra Theorem of SVD ^  88  A2 Detailed Petrographic Descriptions. ^  99  A3 Whole-Rock Chemical Analysis ^  104  A4 Electron Microprobe Analysis . ^  105  A5 MASSBAL.0 source code and test data. ^  120  v  LIST OF TABLES  Table 1.1 Summary of mathematical developments and applications of mass balance modeling in literature. ^  5  Table 1.2 Symbol notation for review of literature (Table 1.1) ^ 8 Table 2.1 Symbol notation for development of mass balance algorithm. ^ 13 Table 2.2 Critical R 2 values for acceptance of solutions used in literature. ^ 19 Table 2.3 AX 2 values as a function of confidence level and degrees of freedom ^ 23 Table 3. la Synthetic chemical data set without assigned errors ^ 33 Table 3.1b Results of calculation with MASSBAL.0 on synthetic data sets without simulated errors ^  33  Table 3.2 Random errors and the synthetic data. ^  36  Table 3.3 Results for synthetic data set No. 1. ^  37  Table 3.4 X2 and R 2 fitting results for synthetic data sets involving three phase fractionation. ^  40  Table 3.5 x2 results of synthetic data sets with inaccurate errors. ^ 43 Table 3.6 x 2 results for the examples with incorrect models. ^ 46 Table 4.1 Chemical data for assimilation problem. ^  52  Table 4.2a Comparison of solutions from MASSBAL.0 and the literature. ^ 53 Table 4.2b Confidence regions around the X2 solutions on the assimilation problem. ^  53  Table 4.3a Chemical Data for Diamond Crater's problem. ^ 58 Table 4.3b The X 2 optimal solutions and the confidence regions on around the X 2 optimal solutions for Diamond Crater's problem. ^  58  Table 4.4 Phenocryst, microphenocryst and groundmass mineralogy in Mount Meager rocks ^  62 VI  Table 4.5 Whole rock chemistry of the Mount Meager volcanic rocks. ^ 65 Table 4.6 Estimates of XRF machine precision and sample preparation uncertainty for sample MM-32. ^  68  Table 4.7 Chemistry of minerals from PLA ^  70  Table 4.8 Chemical data for PLA fractionation problems ^  73  Table 4.9 Mass balance calculation results for PLA. ^  75  Table A3.1 XRF operating conditions for the analysis of major and minor elements. ^  106  Table A3.2 XRF operation conditions for the analysis of trace elements. ^ 107 Table A3.3 The effect of Fe2+/Fe3+ ratio on CIPW norms. ^ 108 Table A4.1 Electron microprobe analyses of PLA plagioclases. ^ 109 Table A4.2 Electron microprobe analyses of PLA amphiboles . ^ 114 Table A4.3 Electron microprobe analyses of PLA pyroxenes. ^ 115 Table A4.4 Electron microprobe analyses of PLA biotites ^ 116  VII  LIST OF FIGURES  Figure 2.1  X2 probability density function ^  22  Figure 2.2 Two dimensional solution spaces defined for two levels of confidence ^ 25 Figure 2.3 Comparison of the goodness-of-fit for different models. ^ 27 Figure 2.4 Flowchart of the calculation procedures of the program MASSBAL.C. ^ 28 Figure 3.1 Procedures to generate error-assigned synthetic data. ^ 35 Figure 3.2 Confidence region around the best solution for error-assigned synthetic example 1 ^  38  x2 and R2 optimization . ^ 41 Figure 3.4 The effect of error estimation on the x 2 technique optimization. ^ 44 Figure 3.3 Influence of standard deviation on  Figure 3.5 Confidence regions on solution for model with "extra" phase. ^ 48 Figure 3.6 Comparison between (false) 2-D confidence regions and projections of (true) 3-D confidence regions. ^  49  Figure 4.1 Two dimensional confidence regions around optimal solutions for assimilation problem (late stage). ^  55  Figure 4.2 Three projections of the 3-D confidence regions around optimal solutions for assimilation problem (middle stage). ^  56  Figure 4.3 Comparison between confidence regions for two phase model and the projections on the corresponding plane for the three phase model for the Diamond Crater's differentiation problem ^  59  Figure 4.4 The location of the MM volcanic complex with the Garibaldi belt of southwestern British Columbia ^  60  Figure 4.5 The element ratios P/K to test the comagmatic relation in PLA rocks ^ 71 Figure 4.6 The projections of three dimensional confidence regions around the optimal solution for the PLA early stage of chemical differentiation. ^ 76  V II I  Figure 4.7 The projections of three dimensional confidence regions around the optimal solution for the PLA late stage of chemical differentiation. ^ 77  IX  ACKNOWLEDGMENTS  This thesis is completed under the supervision of Dr. J. K. Russell, to whom I have much appreciation, for he has been continuously helping and discussing with me during the three years program. His carefully reviewing different drafts of this thesis has improved the manuscript dramatically. Especially, his enthusiasm to science has encouraged me through all the hard times in the research.  I also wish to thank Dr. Tom Brown and K. Fletcher for many good suggestions and comments during both the research and writing periods. I am grateful to the helpful discussions with Bruce James, Yao Cui, Tom Clemo, Xiaolin Cheng, Bo Liang and Min Sun. In addition, I wish to express my appreciation to Stanya for her assistance in major and trace element measurements.  NSERC operating Grant A0820 provided funding for this project.  x  CHAPTER 1 INTRODUCTION AND PRINCIPLES OF MASS BALANCE MODELING This thesis includes the mathematical development of chi square minimization fitting technique (x2 technique) in mass balance modeling and test runs of program MASSBAL.0 using synthetic data sets, in addition to natural data sets. This chapter introduces previous approaches of mass balance calculations in geology and the fitting techniques used in literature.  In Chapter 2, a group of linear equations is developed to express mass balance relations in three phase fractionation process with the amounts of crystallizing mineral phases as the unknowns. Both the residual square minimization fitting technique (R 2 technique) and the x2 technique are discussed, including the limitations of the former and the attributes of the latter. Program MASSBAL.0 was developed to perform all the mass balance calculations, based on the Singular Value Decomposition (SVD), in addition to  X 2 and R 2 minimization principles.  This program was written in C using a UNIX operating system. It provides two results: i) the optimal solution for the mass balance equations and ii) the confidence regions around the optimal solution.  Besides the theoretical derivations, Chapter 3 discusses test runs for program MASSBAL.0 that have been performed to illustrate that the  x 2 technique provides more  reliable solutions than does the R 2 technique. Synthetic examples are used to carry out the test runs, since the expected solutions are known for synthetic data sets. Three groups of synthetic data are involved. First, a synthetic data set without assigned errors is used to test the reliability of the program MASSBAL.C. The second group of data sets is synthetic data with assigned errors. This analysis is used to study that the influence of the standard deviation on the x2 and the R 2 optimizations when the data are subject to measurement errors. The last group of synthetic runs involves incorrect hypotheses (adding the unwarranted phases in the  solution and dropping the essential phases from the solution) in order to test the goodness-of-fit of the x 2 technique.  In addition to the synthetic examples, three natural data sets were also applied to demonstrate that the X 2 technique provides the better solutions than does the R 2 technique. These natural data sets are: assimilation problem in Paricutin volcano, Diamond Crater's fractionation problem and the fractionation processes in Plinth assemblage of Mount Meager volcanic complex.  1.1 Objectives  This thesis is about modeling of geochemical data. The process of modeling data usually includes three steps. The first step involves observation of geological phenomenon, recognition and definition of the problem. Secondly, a model is proposed to explain the phenomenon. The model commonly involves a number of parameters that can be related through mathematical equations. Thirdly, the parameters are fitted by finding the best agreement between the observed data and the values predicted by the model (Press et al., 1987). The array of fitted parameters is termed the "best solution" or "optimal solution". The central part of this thesis emphasizes the third step and utilizes the x 2 technique in mass balance modeling of igneous differentiation. The  X2 technique represents a significant  improvement over the traditional R 2 technique.  Mass, balance modeling has been applied to interpretations of geochemical data since the early 20th century (Bowen, 1928). It is still of importance in quantifying a variety of geochemical processes, especially where other approaches (e.g., thermodynamic modeling) are  2  not possible. Thus, developing a robust and sensitive fitting technique for mass balance modeling remains a significant goal. Three criteria are used to judge a fitting technique. First, the technique must be able to identify a "best solution" for the fitted parameters. If the observed data are subject to measurement errors, the measurement errors should be considered during the search for the best solution. Secondly, the technique must provide a way of evaluating the goodness-of-fit against some useful statistical standards (Press et al., 1987). Lastly, it should also provide a way to estimate the ranges of solution that share similar confidence levels. The optimal solution is usually not the only possible answer since we are working on complicated geological systems. The R2 technique commonly used to search for the optimal solution does not satisfy all three criteria listed above. The specific limitations to the R 2 technique are introduced in the next section and discussed in detail in Chapter 2. The x2 technique developed in this thesis can overcome the deficiencies of the R 2 technique and satisfies the three essential attributes of a mathematical fitting technique.  1.2 Previous Mass Balance Modeling Approaches  Mass balance has been used to solve different geological problems. Metamorphic petrologists employ mass balance and mass transfer concepts to the study of element mobility. They are interested in quantifying the gains and losses of mobile elements during metasomatism and other alteration processes (Myers et al., 1986). Igneous petrologists also use this universal concept to interpret the major element chemical patterns resulting from magmatic differentiation processes, such as fractionation (Nesbitt & Cramer, 1981; Gerlach & Luth, 1980), magma mixing (Wright & Doherty, 1970), partial melting and anatexis (Bea,  3  1989), and assimilation (DePaolo, 1981). Typical problems involve well-defined compositions of the magma end members that can be related by a single process. The mass balance modeling is a reasonable method by which the phases participating in the differentiation process are identified and their amounts are calculated. In the last three decades, there have been approximately 17 papers that have presented either new applications or new mathematical developments of mass balance modeling. The summary of applications, mass balance formulas and fitting techniques in these papers are presented in Table 1.1 with symbol notation in Table 1.2.  Previous applications of mass balance modeling in petrology can be divided into two main groups. The first group of applications includes studies of open geological system processes, such as metasomatic alteration (e.g., Gresens, 1966). The other is the study of closed system processes such as igneous differentiation within a single magmatic system. By closed we mean that material is conserved within the considered object. Fractionation, for example, can be treated as closed system differentiation as long as the compositions of two magma end members and involved mineral phases are known (Stormer & Nicholls, 1977; Wright & Doherty, 1970).  Whereas it is straightforward to apply mass balance to closed system processes, open systems present more difficulty. Where the system is open, we expect mass transfer of at least the mobile elements into or out of the system. However, even in open systems some elements can be recognized as immobile or conserved and separated from those elements which are subject to mass transfer. We consider the systems to be closed for these immobile elements. Mass balance relations of these immobile elements can be used as constraints to calculate the mass gains or losses of the active elements. Therefore, mass balance modeling can still be applied to study open geological system processes. Of course, if no immobile element can be found, the trace of the open system process is very hard to track. 4  Table 1.1 Summary of mathematical developments and applications of mass balance modeling in literature. Authors  ^  Applications  ^  Formulas  ^  Fitting Techniques  Gresens (1966)  Estimate element gains and losses in metasomatic alterations.  PD d i - p i eX i = 100 [ f,(—)  Estimates are based on volume ratios calculated with one immobile element, no 4nalytical errors are considered.  Chakraborty (1977)  Estimate element gains and losses in core-mantle-matrix type metasomatic differentiation.  P•p i =D c •dci + D m •d n: P = Dc + Dm (Mp i =100(D c •dci + D m •dmi - P•p i )/P•p i  Estimates are based on mas s ratio of core and mantle calculated from a pair f immobile elements. Mass ratios from different pa rs of immobile elements are assumed to be consistent, no consideration of analytical or hypothetical errors.  Gerlach & Luth (1982)  Interpretation of differentiation model no formulas for the 1959 Kilauea Iki.  Unweighted least squars fitting.  Francis (1986)  Evaluation of the clinopyroxene paradox.  no formulas  Unweighted least squard fitting, the solution with the R2 less than 0.01 is accepted.  Brophy (1986)  Implications for the origin of Hialumina arc basalt.  no formulas  Unweighted least squart fitting, the critical R 2 value is 0.1.  Grant (1986)  Estimate element gains and losses in metasomatic alteration, another solution for Gresens's equation.  Wright & Doherty Computation for petrologic mixing (1970)^problems.  Pp  eX i = [i; d i - p i ]P  p i = E o i .yi + yn • pi i= 1 n- 1 ,  Graphical solution, cu es of components in original rock versus altered roc are made. Element gains and losses are estimated by omparing curves of mobiles element with curve of immobile element or best fit of a straight line throug a series of points of immobile elements. Unweighted least squarn fitting, the accepted solutions must be non-negative and with the minimum R2.  Table 1.1 continued Authors  ^  Applications  ^  Bartley (1986)  Evaluation of REE mobility in low-grade metabasalts.  Dipple et al. (1990)  Identification of two scales of differential element mobility in ductile fault zone.  Formulas  ^  Fitting Techniques Average volume ratios were estimated by assuming all REEs are immobile, gains and losses of REE were calculated, one whose gain or loss more than analytical error was considered to be mobile element and taken off, new average ratio was calculated again until each REE's gain or loss is with analytical error, they were regarded as immobile. Gains and losses of other mobile REEs were estimated by the average volume ratio of mobile REEs.  eX i = P(fv fpd i - p i )  Immobile elements were used to constrain and calculate the gains and losses of mobile elements.  r PD eX i =100[fv (—d i - p i )] Pp  Unweighted least squares fitting. The solution for which each element satisfies S.<Th • was considered as the best fit.  Bea^Partial melting and anatexis mixing by B=a 0 A.- a i A i (1989)^mass balance among source, i=1 segregate and residue. S 2 = t' I (Si ) 2 i= 1  Stout & Nicholls (1977) Michael & Chase (1987)  Interpretation of fractional crystallization model for petrogenesis of the Snake River Plain Basalts.  Unweighted least squares fitting. n b i .M.^d i E ^ – — (100 +^- p i j^100^100  Studying the influence of primary magma composition, H 2 0 and pressure on Mid-Ocean Ridge basalt differentiation.  no formulas  j=1  Unweighted least squares fitting, acceptable solutions are those with R 2 value less than 0.10 and which makes residuals of each element smaller than analytical uncertainty.  Table 1.1 continued  Authors  ^  Applications  ^  Stormer & Nicholls Interactive testing of magmatic (1978)^differentiation models.  Stout & Nicholls (1983)  Study of origin of the hawaiites.  --'Albarede & Provost Determinations of mode of lunar basalt  Formulas  ^  Mj (b ij -d i )/100^d i -p i  i= 1  no formula.  D.. 0 . Lu•  (1977)  and stoichiometric coefficients for coronitic mineral reactions.  Myers et al. (1986)  Determination of quartz eclogite as a source mixing with carbonate and pelagic sediment as parental Aleutian magmas.  p i = X.d ei^+(1-X-Y)•dsi  Myers et al. (1986/87)  Determination of mixing proportions of two and three components mixing problems with end member compositional variability.  di =  Fitting Techniques Unweighted least squares fitting. Solutions providing the R 2 less than 2.0 and with residual on each oxide less than analytical error are considered better. This method is regarded to work better with the number of components at least as twice big as the number of phases. Unweighted least squared fitting. The solution which has residual squares less than 2.0 is not rejected; no conclusion can be reached for the solution with residual between 2.0 and 3.0 ; The solution with the residual squares more than 3.0 is rejected. Chi-square minimization fitting.  j=1  d.  Graphic solution, C pm -C i curves with different X were made. The X value which can satisfy the variability of all major, trace and rare element was considered as the best fit.  Same as above.  2^I  dI d  Table 1.2 Symbol notation for the summary of literature in Table 1.1. Symbols^Notations P^Mass of parental or original rock D^Mass of daughter or altered rock D m^Mass of mantle rock  D o^Mass of matrix rock M.^Masses of mineral phases Pi^Weight fraction of the ith oxide component in parental or original rock d i^Weight fraction of the ith oxide component in daughter or altered rock b ij^Weight fraction of the ith oxide in the jth phase y i^Weight fraction of the jth mineral in source rock eX i^Element gains or losses from 100 or P grams parental or original rock d ^Weight fraction of the ith oxide component hybrid magma dm^Weight fraction of the ith oxide component mantle rock Weight fraction of the ith oxide component matrix rock de^Weight fraction of the ith oxide component quartz eclogite dei^Weight fraction of the ith oxide component carbonate d:^Weight fraction of the ith oxide component sediment d.^Weight fraction of the ith oxide in mixing component 1 d  2^Weight fraction of the ith oxide in mixing component 2  Density of daughter or altered rock of parental or orginal rock fP^Density ratio between altered and original rocks fv^Volume ratio between altered and original rocks X^Mixing proportion of quartz eclogite or parental magma Y^Mixing proportion of carbonate or component 1 The ratio of the i th element in between the source rock and the liquid. D id^The partition coefficient of the i th element in the j the phase. (Di^The i th phase mass fraction in the whole rock. A^ o^Compositional vector of source rock Compositional vector of segregate B A.^Compositional vector of source minerals (j=1,2,...,n) ^Residual square of oxide component i Si ^ S ^Sum of residual squares Magnitude of imbalance MI PD^  Pp^Density  1.3 Fitting Techniques Used in Literature  A variety of fitting techniques have been applied to determine the parameters in mass balance modeling. These techniques include conserved element fitting (Bartley, 1986), fitting by minimization of sums of residual squares (Stormer & Nicholls, 1977), and fitting by minimization of chi square (Albarede & Provost, 1977).  In order to calculate the element gains and losses between the altered rock and original rock, the volume or mass ratio must be known. Mass balance equations of immobile elements were applied to determine the volume (Bartley, 1986) or mass (Grant, 1986) ratio. This process, which includes: finding immobile elements, estimating the volume or mass ratio between, before, and after the alteration, and calculating element gains and losses, is defined to as conserved element fitting.  Although the conserved element fitting technique is simple to use, it can be biased. There is less problem if there is only one immobile element during the alteration process. When more than one immobile element exists, it has been assumed that there is a single volume or mass change ratio. Generally, one immobile element is used to calculate this volume or mass change ratio (Hanor & Duchac, 1990). Actually, the volume ratios calculated from different immobile elements can be very different due to analytical and hypothetical errors. To avoid bias, every immobile element must be considered during the volume ratio estimation. Actually, least square fitting could probably be used to estimate this ratio using all the immobile elements simultaneously.  Igneous petrologists are more interested in closed system reactions. The parameters of interest are different for different magmatic processes, but commonly include the amounts of 9  mineral phases involved in the fractionation process (Stout & Nicholls, 1977) or the proportions of end members for assimilation (DePaolo, 1981) and magma mixing (Wright & Doherty, 1970) processes. Assimilation and magma mixing can be treated as closed system process, because the end-member chemical compositions are known, as is the composition of the added material.  Mass balance equations for each oxide component are used as constraints to calculate these parameters. Usually there are at least nine such mass balance equations. In most igneous applications the number of unknown parameters is less than the number of equations. Thus, mass balance equations for differentiation processes define an overdetermined linear system of equations. Commonly, the R 2 technique has been used to find the "best solution" (e.g., Stormer & Nicholls, 1977).  As compared with arbitrarily choosing a certain number of equations to estimate the parameters (Olsen, 1984), the R 2 technique is a major improvement, however, it still has limitations. First, it does not consider analytical errors. Secondly, when a solution is achieved by minimization R 2 , there is no criterion for establishing the goodness-of-fit. Last, it does not provide a way to estimate the confidence region around the optimal solution. Thus, a new fitting technique needs to be developed.  The X2 technique developed in this thesis can be dated back to Albarede & Provost's (1977) first study. However, in their work, they mostly emphasize on theory and did not consider the confidence region around the optimal solution. In this thesis, the applications of X 2 technique have been emphasized. The confidence regions around the optimal solutions have  been applied to test different models for three natural data sets.  10  CHAPTER 2 MASS BALANCE MODELING ALGORITHM  This chapter introduces two fitting techniques for mass balance modeling of chemical data. The R 2 technique is reviewed and the x 2 technique is developed. The X 2 is defined as the R2 weighted by the variance (Press et al., 1987). The optimizations of these two techniques are similar where the variances for each oxide are equal. The advantages of the X2 technique over the R2 technique lie in its consideration of analytical errors. The chi-square probability function follows the incomplete gamma function Q(x 2 1 v), which makes it possible to estimate the goodness-of-fit of the  x2 optimization.  The optimal solution to the x2 minimization problem is solved by the Singular Value Decomposition (SVD). SVD is a powerful technique for solving a linear problem with singular or very close to singular matrices. The standard deviations on fit parameters are mutual independent and can be expressed by the confidence regions around the optimal solution that share similar significance (Press et al., 1987). SVD provides a simple way to calculate the confidence regions.  2.1 Mass Balance Relations in Igneous Differentiation  Suppose that two igneous rocks are related through magmatic differentiation. The parent magma is represented by one rock composition p i and the daughter system is represented by the other rock composition d i . The compositions of separating phases are also assumed to be known. In fractional crystallization, for example, compositions of the involved mineral phases may be provided by the analysis of phenocrysts of these two rocks. The following equation expresses the mass balance relation between the parent and the daughter magmas. The mass of the parent magma (P) is equal to the mass of the daughter magma (D) plus the mass(es) of the participating phases (Mi ). 11  P = D + mi . j=1  (1)  P, D and M. in equation (1) represent different phases for different processes. Firstly, for fractionation processes, P and D are masses of two cogenetic rocks. Mi are masses of mineral phases responsible for the change in the system's compositions (Stormer & Nicholls, 1977). Secondly, for assimilation processes, P and D are similarly masses of two cogenetic rocks, whereas Mi are masses of wallrock components (DePaolo, 1981). Thirdly, for magma mixing processes, Mi and D represent masses of end-member magmas (Wright & Doherty 1970). Finally, equation (1) can also be used to describe partial melting and anatexis mixing processes. In this case, P, D and Mi indicate source, residual and segregates, respectively (Bea, 1989).  For all the geological problems being studied, geologists are concerned about the amounts of phases (M.) in order to explain the chemical difference between two magmas. Mass balance calculations are a means to estimate these amounts. A three-phase fractionation case is used as an example to illustrate how mass balance calculations help to solve geologic problems. Symbols used in mathematical development are listed in Table 2.1.  2.2 Fractionation of Three Phases  Two rocks are assumed to be related through fractional crystallization; from petrographic examination, mineral phases Mi (j =1,2,...,m) are considered to have been involved. The question is whether all of the mineral phases or any of the combinations of them can explain the compositional differences between the cogenetic rocks. If so, we also want to know how much of each mineral is needed. In order to answer these two questions, the compositions of the parent and the daughter magmas, as well as the compositions of minerals must be known. 12  Table 2.1 Symbol notations used in development of mass balance modelimz alLwrithm.  Symbols^Notations A^ unweighted compositional matrix A'  aij  [— ]9x3, weighted compositional matrix ai  unweighted compositional array B'  bi [—] 97, 1 , weighted compositional array ai  X^[Mi]3xi, mass array of participating mineral phases a ii^mii - d i , chemical difference between the j th mineral and daughter mauma b i^pi - d i , chemical difference between parental and daughter mamas M.^Mass of fractional crystallized mineral phases 8mj^ranges of weight of fractional crystallized mineral phases P Mass of parental magma Mass of daughter magma D Pi^Weight fractions of the i th oxide component in parental mama d i^Weight fractions of the i th oxide component in dauuhter mii^Weight fractions of the i th oxide component in the j th mineral phase x2 Summed chi squared of nine oxide components ^ R 2^Summed residual squared of nine oxide components RR^Relative error 6.^Standard deviation of analytical error for the i th oxide component x- probability function Q(X21 V)^ Q(01 V)^Limiting value of Q(X 2 1 V) function Q(00 I V)^Limiting value of Q(X 2 I V) function i - j, degree of freedom ^ Difference between certain chi square value and the minimum chi squared value AX2^ p^Probability of ex 2 being less than a certain value U & U'^ixj column-orthouonal matrix VT & V' T^Transposes of an jxj orthogonal matrix V and V' W & W'^jxj diagonal matrix with positive or zero element NP (i)^Column vectors of V' Diagonal elements of W'  13  Usually, the compositions of rocks are determined by X-ray fluorescence analysis and are reported as weight percents of nine oxides Si0 2 , Ti0 2 , Al2 0 3 , FeO, MgO, CaO, Na20, K20 and P2 05 . The analytical errors of the XRF results are also estimated during the analysis. The compositions of the minerals Mi are given by electron microprobe. Equation (1) can be written in terms of the individual oxide constituent: P•pi = D•di + il [ Mi•mii ],^ j= 1  (2)  for i=1,2,...,n, where - pi are weight fractions of oxide component i of the parent composition, - di are weight fractions of oxide component i of the daughter composition, and - mii are weight fractions of oxide component i of mineral phase j.  If the mass of parent system P is defined as 1 mass unit, (1) implies:  ti mj . j=1  D = 1-  (3)  Equation (3) together with equation (2) produces:  rii ivy m i; _ d i) = pi - di.  (4)  j= 1  Equation (4) is the final mass balance equation for fractionation. It can be written in an abbreviated form as (e.g., Stout & Nicholls, 1977; Stormer & Nicholls, 1977): ^ i r l (Mi aii ) = bi. j=1 14  (5)  where - aij are the difference of the weight fractions of oxide component i between the mineral phase j, mij, and daughter's, d i , - b i is the weight fraction difference between parent and daughter magmas, - n is the number of oxide components, and - m is the number of involved mineral phases.  Commonly there are at least nine oxide components in the system that provide nine mass balance equations similar to equation (4). The following group of linear equations expresses the three phase fractionation. all^a12 a13'.\ a21 a22 a23 a31 a32 a33 a41 a42 a43 (M 1\ a51 a52 a 53 M2 a61 a62 a63 \M3) a71 a72a73 a81 a82 a 83 91 a92 a93/  =  (6)  In matrix notation: A * X = B,^  (7)  where A =[a ii ] 9x3 , X = [MORI, and B=Ibiii9 x i.  The matrix A and array B comprise the composition vectors of the two magmas (pi and d) and the involved mineral phases (m ij ). B is the measured chemical difference pi - di between the two rock compositions representing the two magmas. The array X is the unknown vector representing the amounts of minerals M j participating in the fractionation. The left15  hand sides of (6) represent the fractionation model subject to adjustable parameters M. The problem is to fit parameters Mi so that the best agreement between the observed data and the model can be obtained. There are a larger number of equations than unknowns in the linear set of equations (6), thus, it becomes an overdetermined linear problem and can not be solved explicitly. The "best solution" is an approximation that can be obtained by some fitting techniques.  2.3 Optimization of Mass Balance Equations  The most common approach that has been used to find the optimal solution is the unweighted R 2 minimization fitting technique (e.g., Stout & Nicholls, 1977). It has been applied to geology for at least two decades (e.g., Wright & Doherty, 1970). The principle of the R 2 technique is to search for the optimal solution which produces the least R 2 value (R2min), where R 2 (sum of squares of residual) is given by: R2  9^3 = E (b i E (aijMi)) 2 -  i =1^j =1  (8)  The R2 technique represents a quantitative tool in the interpretation of geochemical data, however, it has several limitations. One of its limitations is that the R 2 function does not involve any parameters representing analytical uncertainty. Analytical errors associated with XRF analysis can be as high as 2% relative error for major elements and 8% for trace elements. These errors influence the compositional matrix A and array B in equation (6). Thus, they can affect the "optimal solution" found by R 2 optimization.  16  Recognizing the effects of analytical errors on R 2 optimization, previous workers have tried to solve this problem. In Bea's application of mass balance to partial melting and anatexis processes (1989), he invented a parameter Thj which was related to analytical errors. The solutions which could make the residual for each oxide component less than Thj were considered acceptable. In Stormer & Nicholls's (1977) early mass balance modeling, they simply replaced Thj by analytical errors in their mass balance modeling. They considered the solutions which make the residual for each oxide component less than analytical error were better.  All these previous trials were good. However, they did not solve the fundamental problem; the measurement errors should be considered during the optimization procedure. x2 minimization techniques are designed to handle variances associated with measurement error and are therefore an improvement over the traditional R 2 minimization techniques.  The x2 function can be expressed as (9). The principle of the )( 2 technique is to search the optimal solution which can provide the least )( 2 value (x 2 iiii.).  7C2 =^[ (bi -^(aijMj))2 / 05 i2^(9) i = 1^j=1  where  - )( 2 is the sum of the chi squares of each equations, and - ai is the standard deviation of measurement error of pi.  The )(2 technique is related to the R 2 technique. Its formula is R 2 weighted by the variance ai l. The x 2 technique actually deals with a new linear combination (10):  A' • X ---- B',^  17  (10)  where  ,  r aij^  r bij  = 1. 01 19x3, X = [Mil3x1, and B . =L — i9xi. ,  One of the advantages the X 2 technique has over R 2 technique is that the solution found by X 2 is more sensitive to the constraints imposed by the precise data than to imprecise data. The x 2 technique provides a way to treat good and bad analyses differently. The bigger the analytical uncertainty associated with the oxide component measurement, the smaller its contribution to the x 2 function (9). Thus, the optimal solution from X 2 technique appears to be more reliable solution than the R 2 optimal solution.  In order to get a good solution from the X 2 technique fitting, analytical error must be measured precisely. Both overestimation and underestimation for analytical errors influence the X 2 optimization adversely. The overestimation of error reduces the effectiveness of the data and will make it very difficult to recognize bad hypotheses. Conversely, the underestimation of error gives the data more weight than it deserves and will cause X 2 solutions to tend towards the unweighted R 2 solution. Both poor estimates of error will produce poor solutions.  Not only analytical uncertainty influences the optimal solution, but also errors associated with hypothesis. For example, the number of mineral phases used as the fractionation assemblage can cause a problem. If we keep adding participating mineral phases in the equation (8), we can always gain a smaller and smaller R 2 value. When should we stop adding mineral phases and accept the solution? Different geologist set their own criteria to do this. Their philosophy is that only those solutions which produce a R 2 value less than some critical R 2 value can be accepted. Table 2.2 gives the critical R 2 values used by different  18  Table 2.2 Critical R 2 values for acceptance of solutions.  Authors Francis (1986) Brophy (1986) Michael & chase (1987) Stormer &Nicholls (1978) Olsen (1984) Stout & Nicholls (1977) Bea (1989)  R2 Critical values 0.01 0.1 0.1 2.0 2.0 2.0  Thi *  * Thi is a parameter which is subject to analytical error (Bea 1989).  19  geologists. Compared with these arbitrarily set values, the X 2 technique has another advantage discussed below.  The chi-square probability function P(x 2 1v) is defined as the probability that the observed chi-square for a correct model should be less than a value X 2 (Press et al., 1987). Q(x 2 I v) is its complement and currently used in this thesis. This theorem makes it possible to evaluate the validity of the X 2 optimal solution, including the calculation of the confidence regions around the optimal solution and the goodness-of-fit (Press et al., 1987).  2.4 Mapped Confidence Regions for x2 Optimization Igneous systems are very complicated and errors inevitably exist in the analytical data we use. Therefore, the solution selected by any optimization techniques should not be regarded as the sole answer for the geological phenomenon we observe. There can be many other solutions which are reasonable estimates of the true solution. We use the confidence region around the best solution to represent all the solutions with the same confidence level. Regions defined by the "critical" confidence level should be considered to contain equally valid solutions.  Optimization provides the minimum value of X 2 at a particular set of fit parameters. Obviously, the solutions around this point (x 2 r i n) has larger x 2 values. Let Ax 2 represent the difference between x 2 min and the x 2 value of the observed solution. The 0x 2 value will increase with the observed solution perturbing further away from the optimal solution (Press et al., 1987). The confidence region is the zone around the optimal solution and with a constant 0x 2 value as the boundary (Press et al., 1987). The larger the 0x 2 value is, the larger the confidence region, and the larger the acceptable solution space (e.g. amounts of mineral phases) becomes.  20  Two factors comprise the confidence region. One is the geometric shape, for example, one dimensional problems involving a single fit parameter have confidence regions which are line fragments. Two dimensional problems involving two fit parameters have confidence regions which are ellipses. Three dimensional problems involving three fit parameters have confidence regions which are ellipsoids. Higher dimensional problems involving higher numbers of fit parameters have confidence regions which are polygonal geometry. The confidence regions are centered on the optimal solution. The second characteristic of the confidence region is the size which reflects the confidence level corresponding to a constant value of Ax 2 . Three methods can be used to obtain the Ax 2 value corresponding to a desired confidence level. These methods are:  1) Estimating the delta-chi-square probability density function curves, which are expressed in Figure 2.1; 2) Looking up the Ax 2 distribution tables, which is given as an example in Table 2.3; 3) Computing by program GAMMQ.0 (Press et al., 1987). The source code of program  GAMMQ.0 is included in the disk attached at the back of this thesis. This method is introduced below.  The probability distribution of the variable Ax 2 follows the incomplete gamma distribution (Bury, 1975). Gammq function Q(Ax 2 1v) is defined as the probability that the observed Ax 2 will exceed the Ax 2 for a correct model (Press et al., 1987).  gammq(v/2, Ax 2 /2) = 1-p,^  where - p is the probability or the desired confidence level. 21  (12)  0^4  12 16 20 24 28 32 36 A X2  Figure 2.1 The Chi square probability density functions, for several degrees of freedom.  22  Table 2.3 dx2 values as a function of confidence level and de2rees of freedom.  p 68.3% 90% 95.4% 99% 99.73% 99.99%  1  2  V 3  1.00 2.71 4.00 6.63 9.00 15.1  2.30 4.61 6.17 9.21 11.8 18.4  3.53 6.25 8.02 11.3 14.2 21.1  23  4  5  4.72 7.78 9.70 13.3 16.3 23.5  5.89 9.24 11.3 15.1 18.2 25.7  6 7.04 10.6 12.8 16.8 20.1 27.8  This gamma function (14) has limiting values:  Q(0 1 v) = 1,^Q(...1v) = 0,^(13)  where - v is the degrees of freedom, which is equal to the number of equations (one equation for each oxide) minus the number of fit parameters (mineral phases in fractionation modeling), v =n-m.  The critical confidence level used in this thesis is 95 %. With the 0x 2 estimate as the boundary, the size of the confidence region can be determined. SVD provides an easier way to estimate the confidence region around the optimal solution. The linear algebra theorem, that SVD is based on, is presented in Appendix Al. For the three phase fractionation case, the size of ellipsoidal confidence region is related to the parameters of A's decomposed matrices W and V. The boundary of the ellipsoid is given by:  k 2 = w112(v(i)  • 8N4 1 )2 + w222(v(2) • 6M2) 2 + w332(v(3) • om3)2^(14)  where w 11 , w 22 and w 33 are the diagonal elements of W, V (1) ,V (2) and V(3) are the column vectors of V, 81\4 1 , 61■42 and 3M 3 are the ranges of the weight of three mineral phases.  While Ax 2 is the boundary of the confidence region, the vectors V (i) are unit vectors along the principal axes of the confidence region ellipsoid. The semi-axes have lengths equal to I Ax I /w i • The confidence region for the two dimensional case is illustrated by Figure 2.2.  24  Figure 2.2 Sketches of two dimensional solution spaces defined by two different levels of confidence. The regions are centered on the optimal solution (double circle) and the larger region is associated with higher confidence. Major and minor axes are parallel to unit vectors V(i) and V(2) respectively, and have magnitudes of a/w i and a/w 2 for larger ellipse, b/w i and b/w 2 for smaller ellipse.  25  2.5 Goodness-of-Fit for Optimal Solutions  In the study of magmatic differentiation processes, different hypotheses can be used to explain the geological phenomena. To study fractionation, for instance, we can form either two-phase or three-phase fractional crystallization models to explain chemical differences between parent and daughter magmas. The confidence region can be used to analyze the goodness-of-fit of the x 2 optimal solution, and further test the different models. Optimal solutions with more compact 95 % confidence regions are regarded as better solutions. Figure 2.3 shows that two different models provide different confidence regions. Model b's 95% confidence region is more compact, thus, the optimal solution of model b is considered better than that of model a (Press et al., 1987).  The zero line of one mineral phase is defined as the line on which the amount of this mineral phase involved in is zero (Figure 2.3). If the 95% confidence region includes the zero line of one mineral phase, the solution with zero amount of this mineral phase, is included in the 95% confidence region, and it is possible that this mineral phase did not participate in the magmatic differentiation process.  2.6 Program MASSBAL.0 and Calculation Procedures  Program MASSBAL.0 was written to search for the fitted parameters of igneous differentiation models which purport to explain the chemical differences between two rocks. For example, it calculates the amounts of fractional mineral phases. It does two major calculations. First, the best solutions are calculated using both X 2 and R 2 techniques. Secondly, it calculates the confidence regions around the optimal solutions. The calculations  26  Figure 2.3 The different solution spaces used to compare the goodness-of-fit for different models. The double circles represent the optimal solutions and are centered to the confidence region ellipses. The dashed lines are zero lines (see text).  27  are based on the compositions of the two magma end members and the mineral phases and the procedures can be described as the following steps:  1. Calculate the compositional matrix A =[a ii ] 9x3 , and the compositional array B=[bi]9xi;  2. Convert the relative errors into the absolute values, ai = RR x pi; if RR=O, 6i=1; aij  3. Calculate the weighted compositional matrix A' =f --1 9x3 , and the weighted di  compositional array B''[  bi  ]9x1;  4. Decompose the weighted compositional matrix A' by Singular Value Decomposition (SVD), A' =1.P.WWIT; 5. Calculate the fitted parameters by the X 2 technique; 6. Calculate the minimum X 2 value (x 2inin ); 7. Calculate the significance level of the X 2 optimization through gammp function; 8. Decompose the compositional matrix A by SVD, A =U-W•VT; 9. Optimize the fitted parameters by the R 2 technique; 10. Calculate the minimum R 2 value (R 21nin ); 11. Calculate 0x 2 value which is corresponding to certain probability p, such as 95% by gammp function; 12. Calculate the orientations and lengths of the confidence region.  A flowchart in Figure 2.4 is used to illustrate the calculation procedures of the program MASSBAL.C.  28  compositions of magma end members1  compositions of mineral phases  XRF analytical errors matrix A — — array B SVD weighted matrix A'  weighted array B' X  U, V, W gammq function  2  II^optimal solutionsi  confidence region around optimal solution  Figure 2.4 Flowchart of the calculation procedures of the program MASSBAL.C.  29  This program was written in UNIX based C. The SVD and x2 routines are from Press et al. (1987). The source code, the format of input and output are given in the Appendix. The executable code is included in the attached disk at the end of the thesis.  30  CHAPTER 3 SIMULATED DATA SETS AND ALGORITHM TESTING  This chapter demonstrates how to apply the program MASSBAL.0 to the study of igneous differentiation and compares solutions derived by x 2 and R 2 optimization using synthetic data sets.  3.1 Calculation of Simulated Data Set.  The synthetic data set comprises four synthetic basalts P, D1, D2 and D3, and three standard minerals An 70 , Fo 80 and Di 80 (Table 3. la). The composition of basalt P was synthesized by modifying a natural basalt. The compositions of An 70 , Fo 80 and Di 80 were calculated from their chemical formulas Ca 0.7 Na0.3 A1 1.7 Si2.3 0 8 , Mg 1.6 Fe0.4 SiO4 and CaMg 0.8 Fe0.2 Si 2 06 respectively. The derivative basalt compositions are created by simulating crystallization between a parent (P) and the daughter (D) compositions. The daughter magma D 1 represents crystallization of 5g An 70 out of 100g P. Crystallization of 5g of each An 70 and Fo 80 out of 100g P produces D2. Crystallization of 5g of each An 70 , Fo 80 and Di 80 produces D3.  The normalized compositions of magmas D1, D2 and D3 were calculated from the following formulas:  WtDi = ( 100*Wtp - 5*WtAn70 )/95.0  (14)  WtD2 = ( 100*Wtp - 5*Wt An7c, - 5*WtFoso )/90.0  (15)  WtD3 = ( 100*Wtp - 5* WtAn70 - 5* WtFoso - 5*WtDiso )/ 85 . 0  (16)  where Wt represents weight fraction of oxide component. 31  Program MASSBAL.0 was used to calculate the relationships between the parent magma P and each daughter magma and thereby to explore the effectiveness of the x2 optimization technique. Three models were formed to explain the chemical difference between the magmas. These models were the same as the calculation procedure used to create the data sets: one phase fractionation from P to D 1 , two phase fractionation from P to D2 and three phase fractionation from P to D3.  The comparison between the calculated and the expected solutions are presented in Table 3.1b. The former are close to the true solutions, which are 5.0g for any of the participating mineral phases, within computer roundoff errors. The X 2 and R2 optimal solutions are identical, and the minimum values x 2 „,in and R 2min are both zero because there is no error. All the results are consistent with the expectations for error-free data.  3.2 Simulated Data Sets with Associated Uncertainties.  This section investigates the X2 solutions of synthetic data sets that have associated uncertainties. These synthetic examples are created by adding random errors to the synthetic data set of Table 3.1a.  Equation (18) is used to convert the error-free data into synthetic data with error.  Wt' = Wt + fa^  where - Wt' is the unnormalized new weight fraction of oxide component, - Wt is the error-free weight fraction of oxide component, - a is the simulated standard deviation of weight fraction of oxide 32  (17)  Table 3. la Synthetic chemical data set without assigned errors.  Oxide  SiO 2 TiO 2 Al 2 0 3 FeO MgO CaO Na2 O K2 O P2O5  Minerals  Parent  50.6374* 2.6470 13.6785 11.2741 7.3983 11.2453 2.3189 0.5206 0.2799  An7O  Fo80  Di80  Di  50.5439 0.0000 31.6982 0.0000 0.0000 14.3576 3.4003 0.0000 0.0000  39.1914 0.0000 0.0000 18.7454 42.0632 0.0000 0.0000 0.0000 0.0000  53.1679 0.0000 0.0000 0.0000 7.1330 39.6992 0.0000 0.0000 0.0000  50.6423 2.7863 12.7301 11.8675 7.7877 11.0815 2.2620 0.5480 0.2946  Daughters* D2 51.2785 2.9411 13.4373 11.4854 5.8835 11.6971 2.3876 0.5784 0.3110  D3  51.1673 3.1141 14.2278 12.1610 5.8100 10.0500 2.5281 0.6125 0.3293  * See text for detailed calculation procedures for daughter compositions.  Table 3. lb Results of calculation with MASSBAL.0 on synthetic data sets without simulated errors .  Processes*  Calculated results  Expected results  R2 technique  x 2 technique W. of minerals  x l min  W. of minerals  R2min  P to D 1  M1 = 5.00000g  M 1 = 5.00004g  0.00000  M1 = 5.00004g  0.00000  P to D 2  Mi = 5.00000g M2 = 5.00000g  M 1 = 5.00014g M2 = 4.99998g  0.00000  M1 = 5.00014g M 2 = 4.99998g  0.00000  M 1 = 5.00000g M2 = 5.00000g M3 = 5.00000g  M I = 4.99981g M2 = 4.99994g M 3 = 4.99997g  P to D 3  0.00000  33  MI = 4.99981g M 2 = 4.99994g M 3 = 4.99997g  0.00000  component, and - f is a random number valued between -00 and +00.  The standard deviation (a) is calculated from the relative error by (e):  G = eWt^  (18)  To simplify the modification, f was assigned one of three values: -1, +1 and O. A uniform random number generator was used to generate random numbers ranging from 0 to 1 (Press et al., 1987). Random numbers from 0 to 1/3 caused addition of one standard deviation error to the original weight fraction, 1/3 to 2/3 subtraction of one standard deviation, and 2/3 to 0 no change.  The procedures to generate the synthetic data are illustrated as a flowchart in Figure 3.1. Table 3.2 is an synthetic example with assigned 1% relative errors.  3.2.1 Solutions for Synthetic Data Sets with 1% Relative Error.  The calculated results of the synthetic example from Table 3.2 are presented in Table 3.3. It includes two parts, one part includes the optimal solutions for three fractionation processes, the other part describes the confidence regions around the optimal solutions. The best solutions obtained by the X 2 and R 2 techniques are different and they are no longer identical to the true solutions. The X 2 technique provides a better optimal solution closer than does the R 2 technique. The confidence regions for the two phase fractionation process are illustrated in Figure 3.2. The true solution is included within the 95% confidence region.  34  perfect data Wt  uniform random number generator 0 to 1  Wt x 1 %  0- 1/3 1/3-2/3 2/3-1  f=0 data with error Wt' = Wt + f  Figure 3.1 Procedures to generate the synthetic data sets with assigned errors.  35  Table 3.2a Random numbers to generate error-assigned synthetic data sets. Oxides  Parent P  SiO,  -  TiO 2  +a  A1,0 3  -  a a  Daughters D1 a +a 0 +a -  D3  0  0  0  +a  +a  0 +a -a 0  FeO MgO CaO Na2 O  -a -a  -a  0  -a  K2O P2 O5  0  +a  -a  +a  -a  +a  0  -a  +a  a  -6  0  0  +a  0  -  Table 3.2b Standard deviation calculated from 1% relative error on data in Table 3.1a. Oxide  SiO 2 TiO, A1,0 3 FeO MgO CaO Na.,0 K,0 P,0 5  Parent P  D1  0.5064 0.0265 0.1368 0.1127 0.0740 0.1125 0.0232 0.0052 0.0028  0.5064 0.0279 0.1273 0.1187 0.0779 0.1108 0.0226 0.0055 0.0029  Daughters D, 0.5128 0.0294 0.1344 0.1149 0.0588 0.1170 0.0239 0.0058 0.0031  D3  0.5117 0.0311 0.1423 0.1216 0.0581 0.1005 0.0253 0.0061 0.0033  Table 3.2c Assigned error synthetic data set No. 1. Oxides  SiO, TiO 2 A1,0 3 FeO MgO CaO Na2 0 K,0 P205  Parent P 50.1310 2.6735 13.5417 11.1614 7.3243 11.2453 2.2957 0.5206 0.2771  An 70  Minerals Fo 80  D i 80  D/  50.5439 0.0000 31.6982 0.0000 0.0000 14.3576 3.4003 0.0000 0.0000  39.1914 0.0000 0.0000 18.7454 42.0612 0.0000 0.0000 0.0000 0.0000  53.1679 0.0000 0.0000 0.0000 7.1330 39.6992 0.0000 0.0000 0.0000  50.1359 2.8142 12.7301 11.9862 7.7089 11.1923 2.2358 0.5535 0.2973  36  Daughters D, 51.2785 2.9411 13.5717 11.6003 5.8835 11.6971 2.3876 0.5726 0.3110  D3  51.1673 3.1452 14.2278 12.2826 5.7519 10.0500 2.5028 0.6186 0.3260  Table 3.3 Results of MASSBAL.0 calculations usint! synthetic data set No. I includin i optimal solutions and related confidence bounds.  Processes  Phases  Solutions (g) Expected^x ^R2  R2 min -  -  Confidence regions Semi axes -  Axis orientations  90%^95%^99% P to D i  An 70  5.0  5.459  4.960  9.026  0.081  1.30  1.40  1.59  (-1.00)  P to D,  An 70 Po go  5.0 5.0  4.778 4.850  5.048 5.366  12.566  0.361  1.41 0.66  1.30 0.61  1.62  (-0.992,-0.121)  0.76  (0.121,-0.993)  An70 Po go  5.0 5.0 5.0  5.121 4.966 5.270  5.735 5.565 5.326  5.452  1.45 0.80 0.56  1.58 0.87 0.61  1.82 (0.866,0.156,-0.475) 1.01 (0.497,-0.160,0.853) 0.71 (0.057,-0.975,-0.216)  P to D 3  Di80  37  0.273  Figure 3.2 Confidence regions around the optimal solution for the error-assigned synthetic data set No.1. Bouble circle and cross hair represent the optimal solutions obtained by x2 and R2 techniques, respectively. The open diamond is the true solution and is included within the 95% confidence region.  38  3.2.2 The Effect of a on z2 and R2 Optimization  A group of synthetic data sets with increasing assigned uncertainties were used to compare solutions obtained by the X 2 technique versus the R 2 solutions. The method to create the new synthetic data sets with larger assigned uncertainties is the same as the method described previously except the relative errors used are different.  These simulated data sets are based on a three phase fractionation model which produces D3 from P. Table 3.4 presents the calculated amounts of An 70 , Fo 80 and Di 80 by both fitting  techniques for 11 data sets based on relative errors of 0.1 to 3%. The deviation from the true solution is monitored by the relative deviation (Dev i ), defined by the following equation:  Devi = (^- 5.0 )/5.0,^  (19)  for j = An 70 , Fo g) and Di 80 where - 5.0 represents the expected amount of mineral phases, - M. are the calculated amounts of mineral phases.  Figure 3.3 illustrates the influence of standard deviation on the optimization of the x 2 and the R2 techniques. It shows that as analytical uncertainties increases, the optimal solutions from both fitting techniques deviate from the true solutions. However, the X2 technique continuously presents a much better and more stable solution than the R 2 technique. Within 3% assigned relative error, the optimal solutions from the X 2 technique deviate less than 8% from the expected solution, whereas the solutions obtained with the R 2 technique deviate by up to 42%.  39  Table 3.4 Results for both X 2 and R 2 fittinil of synthetic data sets for three phase fractionation.  £*  Phases  Deviation (%)  Solution (g)  X2  x2^R2  X2min  R2min  R2  0.10  An70 Fo80 Di80  5.015 4.998 5.027  5.074 5.057 5.033  0.306 -0.048 0.542  1.482 1.136 0.658  5.3184  0.0027  0.20  An70 Fo80 Di80  5.025 4.993 5.053  5.148 5.114 5.066  0.496 -0.134 1.060  2.954 2.270 1.310  5.3986  0.0109  0.50  An70 Fo80 Di80  5.065 4.984 5.135  5.370 5.284 5.164  1.298 -0.304 2.696  7.392 5.670 3.274  5.4149  0.0682  0.70  An70 Fo80 Di80  5.089 4.977 5.188  5.517 5.397 5.228  1.782 -0.454 3.768  11.332 7.932 4.566  5.4178  0.1339  1.00  An70 Fo80 Di80  5.121 4.966 5.270  5.735 5.565 5.326  2.428 -0.688 5.600  14.708 11.306 6.504  5.4518  0.2732  1.20  An70 Fo80 Di80  5.146 4.959 5.325  5.881 5.678 5.391  2.920 -0.828 6.502  17.620 13.550 7.420  5.4578  0.3932  1.50  An70 Fo80 Di80  5.186 4.949 5.409  6.098 5.846 5.488  3.622 -1.014 8.182  21.962 16.912 9.758  5.4796  0.6145  1.70  An70 Fo80 Di80  5.211 4.942 5.465  6.242 5.957 5.552  4.216 -1.160 9.292  24.844 19.142 11.044  5.4917  0.7891  2.00  An70 Fo80 Di80  5.247 4.931 5.549  6.457 6.124 5.648  4.940 -1.366 10.988  29.138 22.476 12.968  5.5239  1.0922  2.50  An70 Fo80 Di80  5.310 4.913 5.693  6.812 6.400 5.808  6.150 -1.736 13.852  36.244 28.008 16.160  5.5478  1.7064  3.00  An70 Fo80 Di80  5.370 4.895 5.838  7.164 6.676 5.966  7.002 -2.140 16.756  43.280 33.314 19.328  5.5892  2.4567  * E is the relative error in percent. 40  s.  50  40 0  0  30  a> 7 FE  0  0  0.0^0.5^1.0^1.5^2.0^2.5^3.0^3.5 Model Relative Error (%)  Figure 3.3 Influence of standard deviation a on the optimization of the  41  x2  and R 2 techniques.  The x2 technique is a more robust fitting method than the R 2 technique in light of the effects of random error, however, good estimates of analytical error are essential to proper use of the x2 techniques. Having good analyses implies both precise analyses and good estimation of measurement error. Nevertheless, over- and underestimating the measurement errors inevitably happens in practice. The next exercise investigate the consequences of poor estimation of analytical error on the x 2 optimization. Again the effects are demonstrated through a set of synthetic data which are characterized by inaccurate estimates of error.  The data in Table 3.3 were created with 1% assigned relative errors. When calculating the amounts of An 70 , Fo 80 and Di 8 ,3 needed to produce the magmas D 1 , D2 and D3 from P, instead of using the true relative error of 1%, the calculations assumed other relative errors (e.g. 0.5%). Not all of the relative errors for the oxides were adjusted because the x2 and the R2 techniques are identical when the relative variances of all the oxide components are the same. The results are in Table 3.5 for data sets with assigned relative errors for Si0 2 , Al 2 0 3 , FeO and CaO of between 0.1 and 5%.  Figure 3.4 demonstrates the effects on the optimal solution and x 2 min values with varying relative error. It shows that underestimation of errors causes the solutions to deviate from the true solution much more than does overestimation. Different phases are influenced by the poor estimation of errors differently. However, for An 70 , the solution deviates up to 12% with severe underestimation of error and only 2% for five times overestimation of error. Underestimating errors is the same as ignoring errors in the mass balance calculation, and thus is equivalent to applying the R 2 technique to solve the mass balance problem. On the other hand, overestimation of errors reduces the contribution of certain oxide components to the residuals and thus affects the x 2 optimization. Therefore, the effects of overestimating error are less than the corresponding effects of underestimating error. The mean for a chi-square distribution function is approximately equal to the degrees of freedom of the problem (Bury,  42  Table 3.5 x 2 results of synthetic data sets with inaccurate errors replacirwg the true 1% relative error.  E*(%) Phases  E( % ) Phases  Solution (g)  Dev. (%)  x2min  6.464 6.772 5.547  29.272 35.442 10.938  205.58  1.20  An 70 Fo 80 Di so  5.116 4.956 5.264  -0.882 5.274  5.630 5.650 5.383  12.598 12.996 7.662  98.40  1.30  An 70 Fo - -80 Di x°  5.115 4.954 5.261  2.294 -0.914 5.?17  3.51  5.359 5.286 5.329  7.186 5.720 6.588  49.19  1.50  An 70 Fo so Di so  5.114 4.949 5.255  7.280 -1.018 5.096  2.81  4.976 2.758 6.114  An 70 29.14  1.70  5.115 4.947 5.241  2.294 -1.096 4.824  2.34  Di80  5.249 5.138 5.307  An 70 Pogo Di80  5.194 5.065 5.295  3.888 1.304 5.894  19.26  2.308 -1.080 4.928  2.16  Di so  5.115 4.946 5.246  0.60  An 70 Fo80 Di go  5.164 5.025 5.287  3.282 0.494 5.742  13.72  2.00  An 7o Fo so Di so  5.117 4.945 5.241  2.346 -1.082 4.824  1.88  0.70  An 70 Fo80 Di go  5.146 5.000 5.287  2.916 0.000 5.742  10.33  2.20  An 70 Fo 50 Di so  5.119 4.945 5.236  2.388 -1.108 4.726  1.67  0.80  An 70 Pogo Di go  5.134 4.984 5.277  2.682 -0.318 5.546  8.10  2.30  An 70 Fo 50 Di so  5.120 4.944 5.234  2.410 -1.110 4.678  1.58  0.90  An70 Fo80  5.126 4.973 5.274  2.530 -0.534 5.470  6.56  2.50  A n70 Fo io Di go  5.123 4.944 5.230  2/456 -1.132 4.590  1.44  5.121 4.966 5.270  2.428 -0.688 5.400  5.55  3.00  An 70 F080 Di80  5.128 4.945 5.220  2.566 -1.104 4.396  1.21  5.118 4.960 5.267  2.362 -0.784 5.336  An 70  5.45  5.00  5.143 4.947 5.197  ").86? -1.064 3.936  0.86  0.10  0.20  An 70 Fo80 Di 80 An70 Fo80  Di 80 0.30  An 70 Fog° Di80  0.40  0.50  An 70 Fog°  Din 1.00  1.10  An 70  Fo80 Di go  An 70 F080  Di 80  F080  D i80 An70  1.80  F080  Fo so Diso  * E are the relative errors which replace the true value I %.  43  Solution (g)  Dev. (%) x 2 m i n 2.320  4.00  30  225,  • 25  1751-4  0  20^Pco^125 I- 1  r;  15 0  >.<  75  10  7  2 5fr  5  ^• -25  '  0  1  2^3  Model Relative  4  5^6^1  2^3^4  5  0  Error (7.)^Model Relative Error (%)  Figure 3.4 The effects of poor estimates of standard deviation a on X 2 optimization.  44  1975). The calculated minimum  x 2 value can be compared against the expected mean as a test  of the quality of error estimation. In the problem represented in Figure 3.4 the degrees of freedom are 6. Note that the overestimation of a gives minimum  x2 value of between 0.86 to  4, whereas underestimation of a yields values of between 6 to 205, as the estimates of a get poorer, the x 2 min rapidly deviates from the expected value of F_-- 6.  3.2.3 x2 Confidence Regions: A Measure of Goodness of Fit and A Test of Hypothesis -  -  The error estimation is not the only factor to influence the  X 2 mass balance calculation,  errors in the model or hypothesis also play a major role. In some instances, the number of the participating mineral phases can be uncertain. An extra and unwarranted mineral phase may be added to the postulated fractionating mineral assemblage or, through oversight, an essential phase may not be included in the hypothesis. The calculated confidence region around the optimal solution can be used as a measure of goodness-of-fit and can test the model or the petrologic hypothesis. The effect of error in the model is demonstrated by another set of synthetic test runs.  Two types of synthetic examples based on the models in Table 3. la were made: i) by adding an unrequired phase to and ii) by dropping an essential phase from the mineral crystallizing assemblage. The results are given in Table 3.6. It shows that models that dropped one or more essential phases are associated with large compared to the  x2 values (166 to 1422)  x2 values associated with the true solution (5.45). Additionally, the optimal  solutions deviate significantly from the true solution. On the other hand, the models with added, and unwarranted, extra phases have much smaller minimum x2 values (around 7) and the optimal solutions are very small for the extra added phase.  45  Table 3.6 Results of the x2 technique modeling for the examples with incorrect models.  True Solution x 2 sol.  x z min  Three phase P to D 3  An 70 = 5.121 Fo 80 =4.966 Digo = 5.270  Two phases P to D 3  An70 =8.378  5.45  F080=5.935  313.73  (-0.995^-0.100) (-0.100^-0.995)  1.25 0.61  1.35 0.66  1.55 0.75  Three phase P to D 3  An 70 = 5. 121 Fo go =4.966 Disci = 5 - 27 °  5.45  Two phases P to D 3  An70 =2.256 Di80=7.634  713.60  (-0.894^0.447) (-0.447 -0.894)  1.45 0.84  1.57 0.91  1.79 1.05  Three phase P to D 3  An 70 = 5.121 Fo go =4.966 Diso= 5 - 27 °  Two phases P to D3  Fo 80 =4.334  5.45  D180=7.020  166.48  (0.973^0.230) (-0.230^0.973)  0.60 0.94  0.65 1.02  0.75 1.16  Three phase P to D 3  An 70 = 5.121 Fo 80 = 4.966 Di go =5.270  One phase P to D3  An70 =6.537  1422.48  (1.000)  1.29  1.39  1.59  5.45  One phase P to D I  An 70 =5.459  9.03  Two phases P to D I  An70 =5.577 Fo80=0.219  7.77  (0.981^0.195) (0.195^-0.981)  1.30 0.64  1.41 0.69  1.61 0.79  One phase P to D I  An70 =5.459  9.03  Three phases P to D i  An70 =5.509 Fo80=0.201 Di go =0.130  (0.874^0.192^-0.447) 1.43 (0.450^0.032^0.893) 0.85 (-0.185^0.981^0.059) 0.60  1.56 0.93 0.65  1.80 1.07 0.75  Model  rn  Incorrect Hypothesis Model^x2 sol.  x z min  7.59  Confidence regions Orientations Semi-axes 90%^95% 99%  The zero line is used to test the unwarranted mineral phase. Basalt D 1 is produced for P by fractionation of An 70 . If we add the unwarranted phase Fo 80 to the mineral assemblage, the confidence regions in Figure 3.5 shows that Fo 80 is not needed, because the Fo 80 zero line is included within the 95% confidence region. If all solutions within the 95% confidence region are considered acceptable, then Fo 80 is not needed to produce D 1 from P.  Compactness of confidence regions can be used to evaluate the models which have dropped an essential phase. For example. basalt D3 is produced from P by fractionation of three phases An 70 , Fo 80 and Di 80 . Three two phases models including An70 and Fo 80 , An70 and Di 80 , and Fo 80 and Di 80 were used to reproduce D3 from P. Figure 3.6 is the comparison of each these two phase models' confidence region with the corresponding two phase projection from the three phase model's confidence region. It shows that the projections are more compact than the two phase confidence region. Consequently, the three phase model is better than two phase models.  47  Figure 3.5 Confidence regions around the X 2 optimal solution (double circle) for a (false) model with an extra phase (Fo 80). Cross hair represents the optimal solutions obtained with R2 technique. The open diamond is the true solution. Note that the 95% confidence region includes the Fo 80 zero line.  48  Figure 3.6 The comparisons between the (false) two phase confidence regions (dashed ellipses) and 2-D projections of the (true) three phase (An 70 + Fo 80 + Di 80) confidence regions (solid ellipses) with 90%, 95% and 99% confidence limits, respectively for synthetic fractionation. The double circles are the projections of the optimal solution for the (true) three phase model. The stars are the optimal solutions for the (false) two phase model. 49  CHAPTER 4 NATURAL DATA SETS  Three natural data sets were studied using the X 2 technique. These data sets include one assimilation and two fractionation problems. The Paricutin volcano assimilation problem and Diamond Crater fractionation are somewhat controversial topics and are taken from literature (McBirney et al., 1987 and Russell & Nicholls, 1987).  4.1 Paricutin Volcano: An Example of Assimilation  Paricutin Volcano erupted in Michoacan, Mexico and its chemical differentiation has been studied by different petrologists using different methods. Wilcox (1954), using a graphical method, originally interpreted its chemical variations as resulting from a combination of An 70 and Fo 80 fractional crystallization, and concurrent assimilation of up to 25 weight percent continental crust. Miesch (1979) used a method of vector analysis to examine variations within the same series of analyses and concluded that the compositions evolved through three stages. The first could be explained by fractionation of olivine and plagioclase, the second by fractionation of Ca-poor pyroxene and plagioclase, and the third by assimilation of granite wall rocks. McBirney et al. (1987) used unweighted residual square fitting technique to assess the effects of xenolith assimilation with fractionation of specific combinations of minerals for each of the three major stages. Their calculations suggested olivine fractionation for first stage; loss of orthopyroxene and plagioclase with addition of xenolith for the middle stage; and the late stage is explained by a combination of orthopyroxene crystallization and assimilation of xenolith.  The X 2 technique developed in this thesis is used to study the differentiation processes of Paricutin volcano's middle and late stages chemical variation. It is applied to test the 50  contribution of xenolith to the evolution of Paricutin basaltic andesites. The chemical data and the combined models of fractionation and assimilation are the same as McBirney's (1987) data. The difference between W-47-9 and W-48-5 represents the middle stage differentiation involving plagioclase, orthopyroxene and xenolith. The difference between W-48-5 and W16-5 represents the late stage differentiation involving orthopyroxene and xenolith. Table 4.1 is the chemical data of the assimilation problems in Paricutin Volcano. This paper did not report measurement errors, therefore a 1% relative error and the average XRF error were used. The comparison between McBirney et al.'s (1987) R 2 calculations and MASSBAL.0 x2 calculations is given in Table 4.2a.  Table 4.2a shows that the X 2 optimal solutions are very different from McBirney's (1987) R2 optimal solutions. The assigned analytical errors used also influences the solutions. In the group using an average XRF error, the needed amount of xenolith to contaminate the original magma from the x 2 fitting is more than double that from the R 2 fitting. Additionally plagioclase is predicted to be accumulated rather than fractionated. On the other hand, for 1% relative error, the xenolith needed is a little smaller but still larger than the R 2 solution. The plagioclase is fractionated but the amount is about four times smaller than that predicted by McBirney et al. (1987). The amount of xenolith calculated by MASSBAL.0 is more similar to Wilcox (1954)'s conclusion that 25wt. % xenolith have been added, however, the amounts of olivine and plagioclase are quite different from Wilcox (1954)'s answers (3wt. % 01 and lOwt.% P1).  Generally, the )( 2 fitting suggests that we need a larger amount of xenolith to play a role in the calc-alkaline magma's evolution, and the amounts of plagioclase and orthopyroxene are, relatively, not as important as xenolith.  51  Table 4.1^Chemical data of field assimilation from literature.  Oxides  SiO 2 TiO 2 Al 2 O3 FeO MgO CaO Na 2 O K2 O P205 R2  Daughter W-48-5  Plag  Opx  57.76 0.89 17.18 6.46 5.61 6.90 3.69 1.22 0.29  58.65 0.77 17.42 6.06 4.00 6.41 3.89 1.49 0.30  53.35 0.04 29.51 0.93 0.10 11.80 4.42 0.23  53.86 2.46 13.64 28.53 1.48  Xeno^Calc. daughter  72.01 0.20 14.44 1.78 0.89 2.61 3.96 3.00 0.13  = 0.122  Oxides  Si0 2 TiO 2 Al203 FeO MgO CaO Na2 0 K2 0 P205 R2  Parent W-47-9  Parent W-48-5  Daughter W-16-5  Opx  Xeno  Calc. daughter  58.65 0.77 17.42 6.06 4.00 6.41 3.89 1.49 0.30  59.69 0.80 17.17 5.59 3.71 6.12 3.97 1.66 0.28  53.65  72.73 0.20 14.58 1.80 0.90 2.63 4.00 3.03 0.13  59.69 0.73 17.17 5.66 3.66 6.09 3.90 1.62 0.28  1.75 14.16 27.87 1.56  = 0.022  52  58.67 0.94 17.29 5.92 4.03 6.65 3.93 1.52 0.31  Table 4.2a Comparison of solutions from MASSBAL.0 and the literature. W-47-9 to W-48-5 MASSBAL.0 Solutions  Ave. XRF Errors  Literature  1% Relative Errors  No Errors  1.37g 4.55g -13.02g 196.324  6.37 7.00 -2.76  6.7g 6.6g -8.1g  0.229  0.229  0.122  Plag^-2.78g Opx^3.17g Xeno^-18.65g 2 X min^53.808  R 2 min^0.229  W-48-5 to W-16-5 MASSBAL.0 Solutions 1% Relative Errors  Ave. XRF Errors  Literature No Errors  Opx Xeno x 2min  0.86g -6.23g  0.69g -7.03g  0.21 -8,.33  0.55g -10.67g  91.108  R2min  0.036  118.598 0.036  0.036  0.022  Table 4.2b The confidence regions for the x 2 solutions of the assimilation problem. W-47-9 to W-48-5 Phases  1% relative error Semi-axes Axis orientations 90% 95% 99%  average XRF error Semi-axes Axis orientations 90% 95% 99%  Plag Opx Xeno  0.63 0.68 0.79 1.48 1.61 1.87 1.88 2.04 2.36  0.89 0.96 1.11 (0.265,-0.964,-0.034) 1.46 1.59 1.87 (-0.779,-0.235,0.582) 2.69 2.93 3.38 (0.568,0.127,0.813)  (0.102,-0.986,0.128) (-0.993,-0.094,-0.067) (-0.054,-0.134,-0.989)  W-48-5 to W-16-5 Phases  Opx Xeno  ^  1% relative error^ Average XRF error Semi-axes Axis orientations Semi-axes Axis orientations 90%^95% 99% 90% 95%^99%  0.52 2.20  0.56 2.38  0.64 2.72  (0.994,-0.105) (0.105,0.994)  53  0.78 2.24  0.84 2.43  0.97 2.78  (1.000,-0.016) (-0.016,1.000)  The possible ranges of the amounts of mineral and xenolith are also calculated and presented in Table 4.2b. The confidence regions calculated from 1% relative error are used in the related figures. Figure 4.1 is the 2-D confidence region for the late stage differentiation (Xeno + Opx) in Paricutin volcano. Figure 4.2 is the three projections of the 3-D confidence regions for the middle stage differentiation (Pl+Opx+Xeno) in Paricutin volcano. One of the phenomenon from Figure 4.1 and 4.2 is that the 95% confidence region around the X2 optimal solutions for both the three and two phases assimilation problems do not encompass the solutions given by McBirney et al., 1987 .  4.2 Diamond Craters Volcano: An Example of Fractionation  The Diamond Crater volcanic complex (DCV) is located within the high-alumina basalt petrographic province in Southeastern Oregon and comprises Quaternary alkali olivine basalts (Russell & Nicholls, 1987). The rocks from the Diamond Crater are porphyritic with plagioclase and olivine phenocrysts; plagioclase dominates the phenocryst assemblage (16% 24%). Titanaugite is an important phase but is reported only within the groundmass (Russell & Nicholls, 1987). The chemical compositions of the DCV rocks are very similar, with SiO 2 around 47% weight fraction (Russell & Nicholls, 1987). The slight diversity in chemistry between the DCV rocks is thought to derive from the fractionation of both plagioclase and olivine by Russell and Nicholls (1987). The mass balance relationships are investigated here to further test whether augite could play a role in the chemical differentiation of the Diamond Craters lavas.  All the rock chemistry and mineral chemistry data come from Russell and Nicholls (1987), including the analytical errors for the whole rock wet chemistry. Based on the element ratios (Russell and Nicholls, 1987), the compositions of the two most divergent rocks DCVO45  54  Figure 4.1 Confidence regions around the x2 optimal solution (double circle) for the late stage differentiation of Paricutin volcano. The dashed line is the orthopyroxene zero line. McBirney et al. (1987)'s solution is not included within the solution spaces.  55  4  3  oD  2  1  A. C  o  e  (1.37, 4.55)  McBirney's sol. (6.7, 6.6) -2  2  3^4^5^5^7  PI (g)  -10  O (1.37, -13.02) - 11  -12  o -13  a)  X -t4 - 15  McBirney's sol. (6.6, -16 -2^-1  -  ad)  0^1^2^3  Opx (g) McBirney's sol. (6.7, -8.1)  @ (4.55, -13.02)  Figure 4.2 Three 2-D projections of the 3-D confidence regions (Pl+Xeno+Bi) with 90%, 95% and 99% confidence for the middle stage differentiation of Paricutin volcano. Double circles are the projections of the X 2 optimal solution. McBirney et al. (1987)'s solution is not included within the solution spaces.  56  and DCVO41 were used to represent the compositions for the parent and derivative magmas. The chemical data are reproduced in Table 4.3a.  Three models were formed to investigate the mass balance relation between DCVO45 and DCVO41. The first model is three phase fractionation which include Pl, 01 and Px. The second model is two phase fractionation of P1 and 01. The other model is 01 one phase fractionation model. Table 4.3b gives the calculation results, including the optimal amounts of mineral phases and confidence regions around the optimal solutions.  Table 4.3b shows that the three phase (01+Pl+Px) model has the smallest x 2 values. Figure 4.3 are the comparison between the two phase (01 + Pl) model confidence regions and the two dimension projection (01 + Pl) of three phase model (01 + Pl+ Px) confidence regions. It shows that the semi-axes of the projection from the 95% three phase confidence reigon are 0.29g and 0. 83g for 95% confidence region, respectively, while the semi-axes of the 95% confidence region for the two phase model are 0.28g and 1.27g, respectively, thus, the solution space of the two phase model is more compact. The pyroxene zero lines are not included in the projections on Pl+Px and 01+Px planes from the three phase (01+Pl+Px) confidence regions. Thus, the three phase fractionation model is better and augite is essential to produce the chemical difference within Diamond Crater's lavas.  4.3 Mount Meager: Another Example of Fractionation  4.2.1 Introduction  The Mount Meager (MM) volcanic complex, is located about 150 km north of Vancouver, British Columbia in the Coast Mountains (Figure 4.4), and seated between the Lillooet River and Meager Creek. The Mount Meager volcanic field, together with five other 57  Table 4.3a Chemical data for alkali olivine basalts and minerals at Diamond Craters, Oregon.  Parent DCVO45  Daughter DCVO41  01  P1  Px  Relative error  SiO2 TiO 2 Al203 FeO MnO MgO CaO Na2 O K2 O P2 O5  47.38 1.20 17.30 9.93 0.17 8.81 11.02 2.65 0.27 0.21  47.55 1.46 16.86 10.30 0.17 8.19 10.67 2.94 0.41 0.26  38.65  49.20  17.24 0.27 42.77 0.26  49.64 2.13 2.94 9.76 0.22 13.09 20.78 0.45  0.106 1.667 0.173 2.163 11.765 0.341 0.272 1.132 7.407 9.524  Total  99.21  99.05  99.38  Oxides  30.99 0.48 0.07 15.23 2.94 0.10  99.18  99.01  Table 4.3b The confidence regions for the x 2 solutions of three models for Diamond Crater frationation. These models are three phase model (01 + PI + Py), two phase model (01 + Pl) and one phase model (01). Models  Phases Optimal Sol. (g)  90%^95% 01 P1 Py  4.07 10.98 2.93  01+Pl  01 P1  3.69 8.33  01  01  0.78  01+Pl+Px  Axis orientations  Semi-axes  x2 min  99%  0.24 0.62 1.37  0.26 0.67 1.49  0.29 (0.914 -0.346 0.214) 0.77 (-0.313 -0.263 0.912) 1.70 (-0.259 -0.901 -0.349)  253.09  0.26 1.18  0.28 1.27  0.31 1.45  (0.938 -0.345) (0.345 0.938)  1000.57  0.28  0.31  0.35  (1.000  66.04  58  )  13 m (4.07, 10.98) • (3.69, 8.33)  • •  8 7  1^2^3^4^5^6  01 (g)  7  (10.98, 2.93) 6 5  el  4  oe 3 2  1 0  9^10^11^12 PI (g)  8  7  1  m  13  ^  4  I^..  (4.07, 2.93)  6 5 4 a 3 2  1 1^2^3^4^5^6 01 (g)  Figure 4.3 The comparison between the confidence regions for the two phase (01+Pl) model (dashed ellipses) and the 2-D projections of the three phase (01+Pl+Px) confidence regions (solid ellipses) for the Diamond Crater's fractionation. The confidence regions are with 90%, 95% and 99% confidence limits, respectively. The open circles are the projections of the three phase optimal solutions. The star is the two phase optimal solution. 59  130 0  1 20 °  1250  Queen^ ♦ "O  Sound^  N,%,^  &Cs^e.  \^ \^  52 °  -1-BRITISH COLUMBIA SILVERTHRONE -FRANKLIN GLACIER cetV. bALAL GLACIER v. MEAGER CREEK  • ill  • ELAHO VALLEY  IILMOUNT otc , CAYLEY IPGARICIALDI LAKE *MOUNT \GARIBALDI EXPLORER WATTS POINT PLATE  I^  PACIFIC PLATE  48o  130 0  4, 'I'  •MOUNT BAKER  WASHINGTON GLACIER PEAK  JUAN DE FUCA PLATE 125 °  12ID °  Figure 4.4 The location of the Mount Meager volcanic complex with the Garibaldi belt of southwestern British Columbia (Green et al., 1990).  volcanic fields including Mount Garibaldi, Garibaldi Lake, Watts Points, Mount Cayley and Salal Glacier, forms the north-northwestern extension of Garibaldi volcanic belt. Mount Meager occurs at the northern end of this belt (Green et al., 1988). The Garibaldi belt is the northern continuation of the Cascade volcanic belt and is believed to derive from the subduction of the Juan de Fuca plate beneath the North American plate (Riddihough, 1978). The first study looking at the MM rocks involved the Bridge River pumice and ash (Nasmith et al., 1967). Later, the Meager Creek hot springs were investigated as a potential source of geothermal power, which resulted in detailed mapping (Anderson, 1975 and Read, 1977) of the Mount Meager area. Read (1977) produced a 1:20,000 geologic map of the Meager Creek Geothermal Area. He subdivided the Mount Meager volcanic complex into seven major assemblages based on age, field relations, outcrop and geography (Read, 1977; 1991). In order of decreasing age, they are the Devastator assemblage, Pylon assemblage, Job assemblage, Capricorn assemblage, Plinth assemblage, Mosaic assemblage and Bridge River assemblage. Twenty seven samples were collected by Stasiuk and Russell (1988) which cover these assemblages. Mass balance calculations are applied to investigate the chemical differentiation processes within one of these assemblages, namely, the Plinth assemblage. 4.2.2 Mineralogy and Rock Chemistry Data of the Mount Meager Volcanic Complex  The Mount Meager volcanic rocks include a wide variety of felsic through intermediate to basaltic rocks. However, andesite dominates the volcanic edifice. The mineralogy and petrography of each sample was studied in both hand sample and thin section. All samples are highly porphyritic. The phenocrysts observed in the MM rocks are summarized in Table 4.4. Detailed mineralogy descriptions are given in Appendix A2. 61  Table 4.4 Phenocryst, microphenocryst and groundmass mineralogy in the Mount Meager volcanic rocks*.  Assemblage Name^Phenocrysts^Xenocryst^Microphenocrysts ^ Groundmass crysts  01 Opx Cpx Bi Am P1 Ap^Q^01 Opx Cpx Bi Am P1 Ap Q Z 01 Opx Cpx Bi Am P1 Ap Q Op  ^ Pliocene MM-34 X X X X^x^x^X X X^X^X X ^ x^X X^ x^X X MM-36 x^X X^ X^ Devastator MM-27^X X^ MM-33^X^X^X^X^X^ Pylon^MM-28^X^X MM-35^X^X MM-20^X^X X MM-21^X^X MM-25^X X^X  X^ X^ X^ X^ X^  X X^X X^X X^X X^X X X^X  X X^X  X^X X X^ X^X^X^X X^X^X X X^X^X X X^X X^  X X X Capricorn MM-34^X^X X X^X^X^X X X^ MM-38^X X^ X X^X^X X Plinth^MM-26^X MM-29^X MM-30^X MM-40^X MM-39^X  X^X^X^X X X^ X X X X X X X X X^X^X X X X X^ X^ X^X X X^X^ X^X^X X^X^ X X^X X X^X^X X^X^ X X^X X  Mosaic MM-37 X^X^X^X^X^X^X^X^X^X MM-17 X^X^X^X^X^X^X^X MM-18 X^X^X^X^X^ X^X MM-19 X^X^X^X^X^ X^X  ^  Table 4.4 continued Assemblage Name  ^  Phenocrysts^Xenocryst^Microphenocrysts^Groundmass crysts 01 Opx Cpx Bi Am P1 Ap^Q^01 Opx Cpx Bi Am P1 Ap Q Z 01 Opx Cpx Bi Am P1 Ap Q Op  ^Bridge river MM-kl^X X X^ X^X X MM-16^X^X^ X^X^X MM-2B^X X X^X^X^X X^ MM-32^X X X^ X X^X X X MM-K3^X X X^ X X^X X X MM-7^X X X^ X X^X X MM-2A^X X X^ X X^X X  X^X X  * 01 - Olivine, Opx - Orthopyroxene, Cpx - Clinopyroxene, Bi - Biotite, Am - Amphibole, PI - Plagioclase, Ap - Apatite, Q - Quartz, Z - Zircon, Op Opaque.  Major, minor and trace element concentrations of 27 whole rock samples from MM volcanic complex were made with pressed crushed glass pellets and pressed powder, and measured by X-Ray fluorescence spectrometry. The analyses were made with a Philip PW 1410 X-Ray wavelength-dispersion spectrometer in the Department of Geological Sciences. Major and minor elements include Si0 2 , Al 2 0 3 , Ti0 2 , MgO, Fe 2 0 3 , MnO, CaO, Na 2 0, K2 0 and P 2 05. The analyzed trace elements are Nb, Y, Zr, Sr and Rb. Five duplicate samples were run to estimate the precision of sample preparation. XRF "machine precision" was estimated by running the same pellet five times. The operating conditions for the analysis of major, minor and trace elements are given in Appendix 3.  Ferrous iron was measured volumetrically and ferric iron was calculated from the total iron determined by XRF. Duplicate analyses of ferrous iron were performed to ensure a reproducibility of ± 0.05%. Finally, the absorbed and structurally bound water were measured gravimetrically at 100°C and by the Penfield method, respectively. The operational procedures are given in Appendix 3.  Chemical compositions of the volcanic rocks are listed in Table 4.5. The data include major, minor and trace element concentrations, H 2 0+, H2 0 , loss on ignition (LOI) and -  CIPW normative mineralogy. Estimates of sample preparation and machine precision are reported in Table 4.6.  Much attention has been paid during the XRF analysis in order to reduce the uncertainty related to the machine instability. The counts on the sample used as a monitor were strictly investigated. Only when the counting number of the monitor was within 1% of the mean value were the counts on the three unknowns recorded. Otherwise, this group of sample would be reanalyzed. The results in Table 4.6 also show that uncertainties due to sample preparation are about the same magnitude as the uncertainties due to machine instability, 64  Table 4.5 Whole rock chemistry of Mount Meager volcanic complex.  Sample No. Assemblage* UTM  MM-2A MM-2B BRA BRA 684121 684121  MM-7 BRA 682124  MM-16 BRA 664124  MM-20 PYA 628038  MM-21 PYA 628038  MM-24 PCA 604034  MM-25 PYA 642038  MM-26 PLA 649085  SiO 2 (wt. %) TiO 2 Al 2 0 3 Fe 2 O3 FeO MnO MgO CaO Na2 O K2 O P 2 O5  64.35 0.57 16.57 2.09 2.32 0.08 1.92 4.17 5.10 2.23 0.19  67.69 0.46 16.10 1.06 1.86 0.07 1.36 3.06 5.16 2.41 0.15  68.33 0.47 16.45 1.37 1.68 0.07 1.46 3.33 5.48 2.44 0.15  68.64 0.44 16.18 1.44 1.43 0.07 1.29 3.07 5.32 2.46 0.15  68.68 0.42 16.19 1.42 1.36 0.07 1.25 3.03 5.42 2.50 0.14  67.32 0.52 16.60 1.64 1.62 0.07 1.71 3.67 4.78 2.43 0.18  62.40 0.62 17.41 3.33 1.06 0.07 3.21 5.00 4.73 1.65 0.26  62.76 0.71 17.66 3.09 1.54 0.09 2.67 4.79 5.17 2.07 0.24  67.10 0.54 16.02 1.82 1.40 0.08 1.74 3.28 4.99 2.83 0.23  Total  99.59  99.39  101.22  100.49  100.48  100.55  99.76  100.81  100.03  LOT H 2 0+ H2 O'  1.66 1.34 0.28  2.14 1.96 0.19  0.42 0.31 0.11  0.56 0.40 0.08  0.50 1.04 0.06  1.92 2.15 0.03  2.62 1.26 1.56  1.04 0.59 0.58  0.29 0.23 0.07  Nb (ppm) Y Zr Sr Rb  7 18 180 570 41  7 14 153 449 41  5 15 158 506 50  12 17 157 477 44  14 17 152 0.2 52  23 18 180 595 31  4 14 197 1120 20  26 19 224 769 53  30 19 202 712 5  Q (wt. %) Or P1 Di Hy Mt Hm  13.78 13.16 61.45 2.97 5.28 2.18 0.00 0.79 0.40  18.10 14.24 60.00 0.44 5.15 1.11 0.00 0.64 0.31  16.49 14.13 61.04 1.93 4.07 1.40 0.00 0.64 0.31  18.40 14.37 60.10 0.99 3.73 1.49 0.00 0.61 0.31  17.97 14.59 60.40 1.31 3.39 1.47 0.00 0.58 0.29  18.35 14.23 59.06 0.29 5.28 1.70 0.00 0.72 0.37  13.07 9.70 63.58 1.36 8.14 1.33 1.43 0.86 0.54  10.86 12.03 64.21 2.62 5.94 2.27 0.61 0.97 0.49  17.05 16.62 57.44 1.44 4.33 1.89 0.00 0.75 0.48  Ap  BRA: Bridge River Assemblage; PYA: Pylon Assemblage; PCA: Pliocene Assemblage; CPA: Capricorn Assemblage; DGA: Devastator Assemblage; PLA: Plinth Assemblage; MOA: Mosaic Assemblage.  65  Table 4.5 continued  Sample No. Assemblage* UTM  MM-27 DGA 655079  MM-28 PYA 666082  MM-29 PLA 644089  MM-30 PLA 634109  MM-32 BRA 621119  MM-kl BRA 684121  MM-k3 BRA 684121  MM-33 DGA 621080  MM-34 CPA 610080  Si02 (wt.%) TiO2  P205  60.85 0.85 17.76 2.73 2.77 0.09 2.82 5.08 5.56 1.76 0.36  60.96 0.79 17.79 2.02 3.08 0.10 2.89 5.07 5.08 1.58 0.28  68.74 0.47 15.72 1.65 1.32 0.07 1.42 2.70 5.10 3.06 0.17  60.50 0.84 17.74 3.40 2.02 0.09 2.73 5.06 5.43 1.74 0.36  68.82 0.44 16.09 1.14 1.70 0.07 1.31 3.07 5.30 2.46 0.15  66.90 0.50 16.63 0.95 2.75 0.08 1.59 3.49 5.32 2.22 0.17  68.11 0.46 16.23 0.56 2.63 0.07 1.39 3.17 5.28 2.37 0.15  68.22 0.55 16.29 1.30 1.17 0.04 1.27 1.90 3.02 5.86 0.21  66.64 0.53 16.21 3.72 0.01 0.08 1.94 3.53 4.77 2.48 0.18  Total  100.63  99.62  100.44  99.89  100.55  100.61  100.43  99.82  100.11  0.26 0.40 0.16  0.09 0.28 0.14  0.19 0.28 0.06  0.94 0.90 0.10  0.45 0.28 0.08  1.73 1.48 0.24  0.74 0.75 0.10  3.49 1.90 0.67  0.02 0.10 0.06  Nb (ppm) Y Zr Sr Rb  14 20 231 773 29  4.1 15 177 748 24  12 19 166 533 53  11 18 163 500 58  9.0 16 167 485 41  6 17 170 526 38  4 15 157 508 44  10 26 155 4 116  7 19 163 535 48  Q (wt.%) Or P1 Di Hy Mt Hm II Ap  7.24 10.22 67.06 3.57 7.19 2.81 0.00 1.16 0.74  9.36 9.28 66.34 1.79 9.45 2.10 0.00 1.09 0.58  18.39 17.91 56.24 1.05 3.71 1.71 0.00 0.65 0.35  8.44 10.19 67.07 3.06 5.94 3.12 0.27 1.16 0.75  18.36 14.36 59.74 1.11 4.34 1.18 0.00 0.61 0.31  15.3 12.95 61.91 1.00 6.82 0.98 0.00 0.69 0.35  16.89 13.85 60.35 0.92 6.46 0.58 0.00 0.63 0.31  21.21 34.85 35.40 0.00 3.73 1.37 0.00 0.77 0.44  18.14 14.61 58.17 0.00 5.34 0.00 2.59 0.14 0.38  Al203 Fe2 03 FeO MnO MgO CaO Na 2 0 K2 0  LOI H2 O+  H20  -  66  Table 4.5 continued  Sample No. Assemblage* UTM  MM-35 PYA 608106  MM-36 PCA 592101  MM-38 CPA 572095  MM-39 PLA 638097  MM-40 PLA 639101  MM-18 MOA 604991  MM-17 MOA 642161  MM-19 MOA 603994  MM-37 MOA 582109  5i0 2 (wt.%) TiO 2 Al203 Fe2 O 3 FeO MnO MgO CaO Na2 O K2 O P2 O5  63.15 0.66 16.80 4.15 0.52 0.09 2.37 4.57 5.20 1.87 0.21  65.21 0.63 16.69 3.24 0.88 0.08 2.12 4.26 4.97 2.41 0.27  59.06 0.75 18.88 3.40 2.60 0.12 3.27 5.85 5.15 1.02 0.28  68.10 0.48 16.15 3.34 0.04 0.08 1.68 3.28 5.16 2.55 0.14  66.94 0.52 16.75 1.96 1.60 0.09 1.89 3.62 4.90 2.44 0.17  50.68 1.44 15.56 2.53 7.70 0.16 8.00 8.70 3.97 0.78 0.27  51.33 1.34 18.06 2.47 6.38 0.15 6.06 8.63 4.50 0.74 0.35  49.60 1.49 15.31 2.40 8.10 0.16 8.52 8.55 3.78 0.77 0.28  50.91 1.51 16.57 3.80 5.02 0.14 6.24 9.17 4.42 1.32 0.51  Total  99.59  100.76  100.37  101.01  100.86  99.80  100.01  98.95  99.61  LOT 112 0+ H2 O  0.95 0.83 0.11  1.10 1.10 0.12  1.08 1.01 0.44  0.14 0.08 0.02  0.95 0.98 0.14  0.69 0.92 0.40  0.72 1.58 0.15  0.00 0.53 0.32  0.07 0.34 0.08  Nb (ppm) Y Zr Sr Rb  7 17 202 692 36  5 16 240 882 42  3 14 148 779 15  8 15 138 438 46  6 18 155 483 47  17 16 149 581 17  10 17 196 943 13  19 17 150 537 14  12 18 266 1575 6  Q (wt. %) Ne Or P1 Di Hy Mt Hm Il Ap 01  12.90 0.00 11.03 63.56 3.45 4.81 0.04 2.86 0.92 0.44 0.00  14.88 0.00 14.07 60.02 2.57 4.50 0.91 1.63 0.87 0.56 0.00  7.81 0.00 5.95 70.71 1.43 8.99 3.51 0.00 1.03 0.58 0.00  17.92 0.00 14.85 58.87 0.62 4.27 0.00 2.30 0.18 0.29 0.00  17.18 0.00 14.22 59.70 0.27 5.54 2.02 0.00 0.71 0.35 0.00  0.00 0.00 4.57 57.57 15.14 2.02 2.63 0.00 1.99 0.56 15.52  0.00 0.00 4.32 66.55 10.79 0.04 2.55 0.00 1.85 0.72 13.18  0.00 0.00 4.55 56.52 14.46 1.08 2.51 0.00 2.08 0.59 18.22  0.00 1.87 7.76 57.71 16.46 0.00 3.95 0.00 2.09 1.06 9.10  -  67  Table 4.6a Estimates of XRF machine precision for sample MM-32.  1  2  3  4  5  mean  Std*  E*  67.99 0.43 16.11 2.66  67.96 0.42 16.11 2.67  67.99 0.42 16.13 2.66  68.04 0.42 16.10 2.67  67.99 0.42 16.13 2.67  67.99 0.42 16.12 2.67  0.029 0.0045 0.013 0.0055 ---  0.042 1.06 0.083 0.21 1.0*  0.07 1.30 3.04 5.21 2.40 0.14  0.07 1.29 3.05 5.23 2.40 0.14  0.07 1.27 3.05 5.14 2.40 0.14  0.07 1.28 3.04 5.23 2.40 0.14  0.07 1.29 3.05 5.22 2.40 0.14  0.016 0.0055 0.057 0.0048  1.23 0.18 1.10 0.19  P205  0.07 1.31 3.05 5.30 2.41 0.14  Total  99.38  99.26  99.31  99.23  99.30  Repetition Si0 2 TiO 2 Al 2 0 3 FeO Fe2 03 MnO MgO CaO Na2 0 K2 0  Table 4.6b Estimates of sample preparation uncertainty for sample MM-32.  1  2  3  4  5  mean  Std*  E*  69.07 0.44 16.26 2.99  69.25 0.44 16.19 2.84  67.99 0.43 16.11 2.66  69.10 0.44 16.37 2.78  68.82 0.44 16.09 2.87  68.85 0.44 16.20 2.83  0.50 0.0045 0.11 0.12  o.73 1.02 0.71 4.24 1.0*  P205  0.07 1.32 3.03 5.14 2.49 0.14  0.07 1.27 3.07 5.15 2.45 0.14  0.07 1.31 3.05 5.30 2.41 0.14  0.07 1.32 3.04 5.24 2.44 0.15  0.07 1.31 3.07 5.30 2.46 0.15  0.07 1.31 3.05 5.23 2.45 1.44  0.021 0.018 0.078 0.029 0.0055  1.60 0.59 1.49 1.19 3.80  Total  100.93  100.84  99.38  100.90  100.55  Duplicate pellet Si0 2 TiO 2 Al 2 03 Fe2 0 3 FeO MnO MgO CaO Na2 0 K2 0  * Std = Standard deviation; E = the relative error as a percentage. The relative error of FeO was estimated by the titration analysis.  68  except for SiO 2 , Al 2 0 3 and Fe,2 0 3 , for which the former is about one to 10 times larger than the latter. In mass balance calculations, only the errors brought in during the pellet preparation are used to represent the total analytical errors. 4.2.3 Mass Balance Calculations with the Plinth Assemblage (PLA)  Among the seven assemblages being studied in the Mount Meager volcanic complex, the lavas which comprise the Plinth assemblage (PLA) show potentially comagmatic relationships based on the field relations and mineralogy. Mass balance calculations are applied to further study the relations between the PLA rocks and investigate the differentiation processes. 4.2.3.1 Comagmatic Relationships  Five samples were collected from the Plinth assemblage which are numbered as MM-26, MM-29, MM-30, MM-39 and MM-40. All five rocks are highly porphyritic and with large amounts of phenocrystal plagioclase and quartz. The plagioclase is unzoned. Quartz occurs as broken crystals and may be xenocrystal. Medium-sized biotite exists in all samples and varies from a few grains to 10% in volume. Biotite is breaking down to amphibole, opaque and pyroxene. There are also small amounts of pyroxene microphenocrysts in MM-29, MM-39 and MM-40, and a small amount of amphibole microphenocrysts in the groundmass of MM-26 and MM-29 (Table 4.4). The rock chemistry of PLA rocks is presented in Table 4.7. All the samples have very high SiO2 (-2:-. 68 wt. %), except MM-30 which contains about 60 wt. % SiO2. They are classified as rhyodacites, except for MM-30 which is andesite. The element ratios P/K are close together and within analytical error of each other, except MM-30 (Figure 4.5). This similar chemical signature, together with the similar mineralogy suggests that except MM-30, 69  Table 4.7a The average electron microprobe analyses of PLA biotites and amphibolse.  Sample  MM-26 Bi  MM-29 Bi  MM-30 Bi  MM-39 Bi  MM-40 Bi  MM-26 Am  MM-29 Am  Al203 FeO MgO TiO2 MnO Na2 0 CaO K2 0 F Cl  37.16 13.43 17.37 13.91 0.37 4.27 8.84 0.03 0.56 0.17 0.09  36.97 13.36 16.72 14.19 4.30 0.33 0.55 0.07 8.56 0.28 0.10  36.57 13.21 16.60 13.75 4.21 0.40 0.65 0.14 8.32 0.28 0.10  36.69 13.58 17.05 13.96 4.31 0.35 0.60 0.01 8.88 0.51 0.08  40.83 12.22 13.40 14.03 3.04 0.32 1.73 7.75 3.20 0.16 0.01  48.94 5.93 11.46 16.54 0.97 0.78 1.43 11.28 0.32 0.21 0.07  44.88 9.35 10.40 16.03 1.98 0.62 2.14 11.87 0.46 0.17 0.01  Total  96.20  95.43  94.23  96.02  96.69  97.93  97.91  5  8  4  3  6  1  4  Si0 2  N*  Table 4.7b The average electron microprobe analyses of PLA feldspars and pyroxenes.  MM-26 PI  MM-29 P1  MM-30 P1  MM-39 P1  MM-40 P1  MM-29 Py  MM-39 Py  MM-40 Py  Si0 2 Al2 03 FeO MgO MnO Ti 02 Na2 0 CaO K2 0  62.57 23.46 0.22 0.00  61.44 23.92 0.18 0.01  60.62 24.69 0.21 0.01  61.44 23.92 0.18 0.01  62.52 23.50 0.20 0.01  7.92 5.10 0.71  7.63 5.56 0.67  7.10 6.52 0.60  7.63 5.56 0.67  8.01 4.89 0.77  51.12 2.59 5.68 16.11 0.14 0.53 0.32 22.17  54.28 1.71 13.74 28.71 0.37 0.19 0.04 1.31  53.82 2.04 13.66 28.66 0.38 0.16 0.02 1.24  Total  99.98  99.41  99.75  99.41  99.90  98.66  100.35  99.98  4  12  16  25  12  2  3  3  Sample  N  *N is the number of measured points.  70  0.16  T 7-  MM- 3 0  0.12  "*---  0 08 MM-26 MM-40 0.04 MM- 29  MM-39  I 0.00 ^ ^ 16^18 20^22^24^26 28^30 1^1^  Si/K  Figure 4.5 The element ratios of PLA rocks in Mount Meager volcanic complex. Four samples have the P/K ratios close enough within the analytical errors and are comagmatic.  71  the PLA samples are comagmatic and related by fractionation. Although the presence of quartz xenocrysts might suggest that assimilation has played a role, the lack of other evidence prevents further study of assimilation. The mass balance calculations using x2 fitting techniques is applied to study the chemical relationships between these rocks.  The compositions of the rocks are presented in Table 4.5. The compositions of the phenocrysts appearing in these five rocks were analyzed by electron microprobe. The analyzed minerals include plagioclase and biotite for all samples, pyroxene for MM-29, MM39 and MM-40, and amphibole for MM-26 and MM-29. The analyses, calculated structural formula and mole% end members are presented in Appendix A4. Table 4.7 contains the average compositions of the phenocrysts in these five samples.  4.2.3.2 Mass Balance Calculations  In terms of SiO 2 wt%, the chemical difference between MM-40 and MM-26 is considered as the early stage in the chemical differentiation of PLA group and MM-29 and MM-39 as the late stage. The fractionation of both plagioclase and biotite was considered to play a major role in the evolution of these rocks. The other two phases, amphibole and pyroxene exist in the groundmass. We want to find out whether pyroxene and amphibole are crucial to explain the observed chemical differences. Mass balance calculations are applied to carry out this job. Two models are formed. The first one is to test whether amphibole fractionated at the early stages to produce MM-26 from MM-40. The second task is to test the contributions of both the pyroxene and amphibole to the late stage chemical evolution in PLA group producing MM-29 from MM-39. Table 4.8 gives the chemical compositions of the end member rocks and mineral phases for these two processes.  72  Table 4-8 Chemical data of the fractionation problems in MM modeled by mass balance calculations. a. MM-40 to MM-26 Oxides  Si0 2 (wt.%) TiO 2 M203 Fe2 0 3 FeO MnO MgO CaO Na2 0 K2 0 P205  Parent MM-40 66.94 0.52 16.75 1.96 1.60 0.09 1.89 3.62 4.90 2.44 0.17  Daughter MM-26  Mineral phases plagioclase^biotite  67.10 0.54 16.02 1.82 1.40 0.08 1.74 3.28 4.99 2.83 0.23  62.57 23.46 0.22 0.00 5.10 7.92 0.71  amphibole  37.16 0.37 13.43  48.94 0.97 5.93  17.37 4.27 13.91 0.03 8.84 0.56  11.46 0.78 16.54 11.28 1.43 0.32  b. MM-39 to MM-29 Oxides  Si0 2 (wt.%) TiO 2 Al203 FeO MnO MgO CaO Na2 0 K2 0 P205  Parent MM-39 68.10 0.48 16.15 3.05 0.08 1.68 3.28 5.16 2.55 0.14  Daughter MM-29 68.74 0.47 15.72 2.80 0.07 1.42 2.70 5.10 3.06 0.17  plagioclase 61.44 23.92 0.18 0.01 5.56 7.63 0.67  73  Mineral phases biotite^amphibole 36.97 4.30 13.36 16.72 0.33 14.19 0.07 0.55 8.56  44.88 1.98 9.35 10.40 0.62 16.03 11.87 2.14 0.46  pyroxene 51.12 0.53 2.59 5.68 0.14 16.11 22.17 0.32  Plagioclase, biotite and amphibole were considered to participate in the early stage of fractionation. Table 4.9 is the calculated results. It shows that about 6% plagioclase, 2% amphibole and 1% biotite crystallized from MM-40 to produce MM-26. Figure 4.6 contains the three 2-D projections for the 3-D confidence regions. It shows that the zero line of biotite is included in the 95 % confidence region around the optimal solution on the projections on Pl+Bi and Bi+Am planes. Thus, biotite did not contribute to the early stage chemical differentiation in PLA group. It was replaced by amphibole, which is consistent with the petrographic observation that biotite is breaking down to amphibole.  Plagioclase continues to dominate the fractionation assemblage during the late stage chemical differentiation in PLA group. Pyroxene microphenocrysts start to appear and scattered throughout the groundmass. The fractionation models formed by mineral assemblage of plagioclase, biotite, amphibole and pyroxene are used to explain the chemical difference between MM-39 and MM-29. The calculated results are also presented in Table 4.9. It shows that the four phase model (Pl+Bi+Am+Px) has to be rejected since it produces a large minimum X 2 value (6199). Figure 4.7 shows the projections on Am+Px and Pl+Px planes from the three phase (Pl+Am+Px) confidence regions. The zero lines of pyroxene is included in the projections of 95% confidence regions. Figure 4.7 also shows the comparison of two phase (P1+Am) confidence regions and the projection on Pl+Am plane from the three phase (Pl+Am+Px) confidence region. The semi-axes of the former are about 0.50g and 1.60g, respectively, while the semi-axes of the latter are 0.48 and 1.81g, respectively, therefore, the confidence regions for the two phase model are a little more compact. Thus, pyroxene is not considered as a major mineral phase in the late stage of chemical differentiation in PLA group. Only two phases, which are about 10% plagioclase and 3% amphibole caused the difference between MM-39 and MM-29.  74  Table 4-9 The mass balance results of PLA fractionation. a. Results of three fractionation models for MM-40 to MM-26 (one three phase model, two two phase model).  Models  Phases Optimal Sol. (g)  X 2 min  90%  Semi-axes 95%  99%  Axis rientations  1.86 0.97 0.49  2.00 1.05 0.53  2.28 1.19 0.61  (-0.984 -0.156 0.080) (-0,175 0.834 -0.523) (0.015 -0.529 -0.849)  P1 Bi Am  5.05 -1.10 2.47  Pl+Bi  PI Bi  6.06 0.77  504.46  1.91 0.76  2.05 0.82  2.32 0.92  (-0.996 -0.090) (0.090 -0.996)  Pl+Am  P1 Am  5.58 2.00  345.59  1.87 0.58  2.01 0.63  2.27 0.71  (1.000 -0.007) (-0.007 -1.000)  Pl+Bi +Am  325.54  b. Results of two models for MM-39 to MM-29 (two two phase model, two three phase model and one four phase model). Models  Phases Optimal Sol. (g)  X 2 min  90% PI Bi Am Px  -18.59 -9.97 -1.74 4.98  Pl+Am+Px  P1 Am Px  12.153 4.078 -0.741  Pl+ Am +Bi  PI Am Bi  -19.82 -13.83 7.06  P1 Bi  -0.62 -8.24  P1 Am  10.595 3.178  P1+Bi+Am+Px  P1+ Bi P1+ Am  Semi-axes 95% 99%  Axis orientations  0.03 0.27 2.44 3.60  0.04 0.30 2.66 3.92  0.04 0.34 3.07 4.53  (0.062 -0.833 -0.494 -0.242) (-0.139 0.414 -0.301 -0.847) (0.870 -0.069 0.377 -0.311) (-0.460 -0.360 0.723 -0.356)  51.70  0.25 2.92 0.92  0.27 3.16 1.00  0.31 3.62 1.15  (0.058 0.543 0.838) (0.884 0.362 -0.297) (0.465 -0.758 0.459)  6315.52  0.04 2.97 0.59  0.04 3.21 0.63  0.05 3.68 0.73  (0.064 -0.859 -0.508) (0.961 0.191 -0.201) (0.270 -0.475 0.838)  2.04 0.05  2.20 0.05  2.51 0.06  (-0.997 -0.075) (0.075 -0.997)  1.68 0.46  1.81 0.48  2.06 0.56  (1.000 -0.024) (-0.024 -1.000)  6199.72  7313.89 58.47  75  Figure 4.6 Three 2-D projections of the 3-D confidence regions (Pl+Am+Bi) with 90%, 95 % and 99% confidence for the PLA early stage differentiation of Paricutin volcano. The dashed ellipses are the confidence regions with 90%, 95% and 99 % confidence limits for the two phase (Pl+Am) model. Double circles are the projections of the x 2 optimal solution.  76  (4.08, -0 .74)  • 3  be I (11.. —1  -3 -5^......... -1^1^3^5^7  ,  Am (g) 5  ( 1 2.15. -0.74)  • 3  -3 -5  7  2  •  (10.80, 3.18)  .  (12.15, 4.08) 8^10^12^14^lg PI (g)  Figure 4.7 The 2-D projections (solid ellipses) of 3-D (Pl+Am+Px) confidence regions with 90%, 95% and 99% confidence limits, respectively for PLA late stage differentiation. The dashed ellipses are the two phase (Pl+Am) confidence regions with 90%, 95% and 99% confidence limits, respectively. The horizontal dashed lines are the pyroxene zero line.  77  CHAPTER 5 SUMMARY AND DISCUSSION  A new x 2 minimization fitting technique has been developed in this thesis to study mass balance relations during igneous differentiation processes. The X 2 fitting differs from the conventional R 2 minimization by weighting the residual squared with the variance of the measurement error. When analytical data are error free or the variances for all the data points are the same, these two fitting techniques provide identical results. However, when the chemical data (individual oxides) vary in analytical quality, the X 2 is a better fitting technique than the R 2 technique.  The advantages the X2 fitting has over the R 2 fit can be summarized on three points. First, in the X 2 technique, the fitted parameters depend more on the better quality data than on the poorer quality data , while the R 2 techniques treat each data point equally. Secondly, the  X 2 technique can estimate a confidence region around the optimal solution. And third, it is possible to estimate statistically the goodness-of-fit in the X 2 fitting, whereas in the R 2 technique, a critical R 2 value is arbitrarily set to quantitatively assess the acceptability of the optimal solution can be accepted or not.  It was demonstrated using synthetic data that the X 2 fitting is a more reliable technique than the R 2 technique. The optimal solution found by the X 2 technique is much closer to the true solution than that by the R 2 technique. With the relative error for all the data points increasing up to 3% in the synthetic data, the X 2 technique continuously produces a much better and more stable solution than the R 2 technique. Within the range of 3% relative error, the optimal solutions from the X 2 technique deviate from the true solution less than 8 %, whereas the optimal solutions from the R 2 technique deviate more than 42%. Therefore, as measurement errors get greater, X 2 fitting provides a better solution compared to that chosen by R 2 minimization.  78  The importance of the X2 technique is also reflected in the calculated confidence regions around the optimal solutions. The geological systems we study are highly complicated and the analytical data used to describe the system are subject to measurement errors. There is a random component included in the measurement error, thus the optimal solution is not the unique realization of the true solution. Instead, there are many other solutions which can also represent acceptable estimates of the true solution. The X 2 technique brings in the confidence region to represent the possible realizations to the true solution, rather than using a unique optimal solution as is done in R 2 fitting.  The confidence region is the zone of the possible solutions which share the same confidence limit. It includes two factors. One is the shape of the region, and the other is the size reflecting the confidence level. The first one depends on the dimension of the solved problem, which is the number of fitted parameters. The shape of the confidence region is a line segment for one dimension centered on the optimal solution, an ellipse for two dimensions and an ellipsoid for three dimensions. The confidence levels used in this thesis are 90%, 95% and 99%, which represent the probabilities of the true solution falling within the calculated regions (90, 95 and 99%, respectively).  The critical confidence level we apply is 95%. The solutions within the 95% confidence region around the x2 optimal solution are considered to be good estimates of the true solution. For all the synthetic data sets, the true solution is included in the 95% confidence region of the optimal solution. The goodness-of-fit is also investigated through the 95% confidence region. The model with the optimal solution whose 95 % confidence region is more compact is regarded as better model to explain the differentiation processes. For example, if two phase fractionation model has an optimal solution whose 95% confidence region is more compact  79  than that of the three phase fractionation model, then the two phase model is better to explain the chemical difference of two magma end members.  The x2 minimization fitting appears to be a good technique for quantitative interpretation of chemical data. It is powerful to test hypotheses explaining chemical difference of two magma end members. However, the estimation of analytical uncertainties is very important to get a good solution. The under- or overestimation of errors have a large effect on the optimal solution and the confidence regions. The underestimation of error has larger effect on the solution than the overestimation of error. The overestimation of error would bring us a smaller and unreliable confidence region.  Program MASSBAL.0 was written to do all the x 2 mass balance calculations. It uses Singular Value Decomposition (SVD) to search the optimal solution and estimate the confidence region around the optimal solution.  The newly developed x 2 minimization fitting technique has been applied to study three field problems. One is the assimilation in Paricutin volcano of Mexico, the other two are fractionation in Mount Diamond Crater of Oregon and Mount Meager Plinth Assemblage (PLA) group of British Columbia. The magmatic differentiation models of the first two were investigated to clarify some debate on them. The chemical difference of PLA group was used to provide some insight the magmatic evolution of the Mount Meager volcano.  In the field assimilation problem, the x 2 technique was applied to test the differentiation models explaining the chemical variation in the volcanic rocks of Paricutin volcano in Michoacan, Mexico. The evolution of this basaltic andesite to andesite volcano was divided into three stages. The early stage involved the fractionating of 01 and Pl, the middle stage involved the fractionating of Opx and Pl., and the last stage involved the fractionating of Opx.  80  However, the better explanation for the chemical variation can be obtained if a certain amount of xenolith is added in the magmatic system concurrent fractionation of mineral assemblage. The x2 mass balance calculations are applied to test the influence of xenolith assimilation. The X 2 optimal solution calculated from MASSBAL.0 were different from the McBirney et al.'s (1987) results. The X 2 calculations suggest that for the middle stage a larger amount of xenolith (about 20 wt. %) was needed to contaminate the magma on its way to the surface with about 3% plagioclase added instead of taken away from the magma and 3% orthopyroxene fractionation. In the late stage, about 6% of xenolith was added in the magma while about one percent of orthopyroxene crystallized from the magma. The difference between the x2 calculations and the R 2 calculations is so large that the solution space of the 95% confidence region around the X 2 optimal solution does not include the R 2 optimal solution obtained by McBirney et al. (1987). The x2 optimal solution for the amount of xenolith participating the middle stage is consistent to Wilcox's (1954) explanation, which was about 20% xenolith contaminating the magma.  From the Diamond Craters x 2 calculation results, about 3% weight fraction of augite is essential in addition to 11 % plagioclase and 4% olivine, to produce the slight chemical difference within the DCV basalts. This result is different from Russell and Nicholls (1987) who concluded that only P1 and 01 were needed. The existing of augite in the groundmass actually infers that the augite has already participated in the fractionation process.  The chemical difference of PLA rocks were divided into early and late two stages based on the petrography and element ratios. Three phase plagioclase, amphibole and biotite fractionation model was used to explain the early stage chemical evolution. The zero line of biotite is included within the 95% confidence region around the x2 optimal solution, thus, very possibly biotite has stopped crystallizing to produce MM-26 from MM-40. In the late stage,  81  only plagioclase and amphibole participated the fractionation process to produce MM-29 from MM-39, while pyroxene is not needed. Although we can not prove that the X 2 solution is better than the R 2 solution, the x2 technique provides the confidence region around the best solution. We can estimate the goodness-of-fit by investigating the compatibility of confidence region, which has not been tried before. The further improvement in the x 2 fitting can be obtained by more precise measurement errors and the representative compositions of mineral phases and xenolith.  This thesis is the first try to use X 2 as a maximum likelihood estimate in modeling geological data. It reveals the importance of quality in analytical measurements where those chemical compositions are to be used in modeling. Also it treats differently the good and bad analyses, thus accentuating the need for precise analyses and the reporting of all analytical errors.  82  REFERENCES Albarede, F., and Provost, A., Petrological and geochemical mass-balance equations; an algorithm for least-square fitting and general error analysis, Comput. Geosci., 3, 309-326, 1977. Anderson, R.G., The geology of the volcanics in the Meager Creek map area, southeastern British Columbia, B. Sc. thesis, 130pp., Univ. of British Columbia, 1975. Bartley, J.M., Evaluation of REE mobility in low-grade metabasalts using mass-balance calculations, Norsk Geologisk Tidsskrift, 66, 145-152, 1986. Bea, F., A method for modelling mass balance in partial melting and anatectic leucosome segregation, J. Metam. Geol., 7, 619-628, 1989. Bowen, N.L., The evolution of the igneous rocks, Princeton University Press, 332pp., 1928. Bury, K.V., Statistical models in applied science, New York, Wiley, 625pp, 1975. Chakraborty, K.R., A new method of material balance evaluation in metamorphic differentiated systems, Contrib. Mineral. Petrol., 65, 101-110, 1977. DePaolo, D.J., Trace element and isotopic effects of combined wallrock assimilation and fractional crystallization, Earth and Planetary Science Letters, 53, 189-202, 1981. Dipple, G.M., Wintsch, R.P., Andrews, M.S., Identification of the scales of differential element mobility in a ductile fault zone, J. Metam. Geol. , 8, 645-661, 1990. 83  Francis, D., The pyroxene paradox in MORB glasses -- a signature of picritic parental magmas, Nature, 319, 586-589, 1986. Gerlach, T.M., Luth, W.C., Mass balance differentiation models for the 1959 Kilauea lid lava lake, Eos Trans. AGU, 61, 1143, 1980. Grant, J.A., The isocon diagram; a simple solution to Gresens' equation for metasomatic alteration, Econ. Geol., 81, 1976-1982, 1986. Green, N.L., R.L. Armstrong, J.E. Harakal, J.G. Souther, and P.B. Read, Eruptive history and K-Ar geochronology of the late Cenozoic garibaldi volcanic belt, southwestern British Columbia, Geol. Soc. Am. Bull., 100, 563-579, 1988. Gresens, R.L., Composition - volume relationships of metasomatism, Chem. Geol., 2, 47-65, 1966. Hanor, J.S., and K.C. Duchac, Isovolumetric silicification of early Archean komatiites; geochemical mass balances and constraints on origin, J. Geol., 98, 863-877, 1990. McBirney, A.R., H.P. Taylor, and R.L. Armstrong, Paricutin re-examined: a classic example of crustal assimilation in calc-alkaline magma, Contrib. Mineral. Petrol., 95, 4-20, 1987. Michael, P.J., and R.L. Chase, The influence of primary magma composition, H 2 O and pressure on Mid-Ocean Ridge basalt differentiation, Contrib. Mineral. Petrol., 96, 245-263, 1987.  84  Myers, J.D., and C.L. Angevine, Mass balance calculations with end member compositional variability: application to petrologic problems, Eos Trans. AGU, 67, 404, 1986. Myers, J.D., C.L. Angevine, and C.D. Frost, Mass balance calculations with end member compositional variability: applications to petrologic problems, Earth Planet. Sci. Lett., 81, 212-220, 1987. Myers, J.D., C.D. Frost, and C.L. Angevine, A test of a quartz eclogite source for parental Aleutian magmas: a mass balance approach, J. Geol., 94, 811-828, 1986. Nasmith, H., and G.E.R. Mathews, Bridge River ash and some other recent ash beds in British Columbia, Can. J. Earth Sci., 4, 163-170, 1967. Nesbitt, H.W., and J.J. Cramer, Graphical representation of mineral equilibria and materials balances in igneous rocks, Contrib. Mineral. Petrol., 78, 136-144, 1981. Nicholls, J., The statistics of Pearce element diagrams and the Chayes closure problem, Contrib. Mineral. Petrol., 99, 11-24, 1988.  Olsen, S.N., Mass-balance and mass-transfer in migmatites from the Colorado Front Range, Contrib. Mineral. Petrol., 85, 30-44, 1984.  Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling, Numercal recipes in C, Cambridge University Press, 818pp., 1987.  Read, P.B., Geology of Meager Creek geothermal area, British Columbia; Geological Survey Canada, Open File 603, 1977. 85  Read, P.B., Mount Meager complex, Garibaldi belt, southwestern British Columbia, Geosci. Can., 17, 167-170, 1991.  Riddihough, R.P., The Juan de Fuca plate, Eos Trans. AGU, 59, 836-842, 1978. Russell, J.K., and J. Nicholls, Early crystallization history of alkali basalts, Diamond Craters, Oregon, Geochim. Cosmochim. acta, 51, 143-154, 1987. Stasiuk, M.V., and J.K. Russell, Petrography and chemistry of the Meager Mountain volcanic complex, southwestern British Columbia, Current Research, Part E, Geol. Survey Can., Paper 89-1E: 189-196, 1989.  Stasiuk, M.V., and J.K. Russell, Quaternary volcanic rocks of the Iskut River region, northwestern British Columbia, Current Research, Part E, Geol. Survey Can., Paper 90-1E: 227-233, 1990. Stasiuk, M.V., and J.K. Russell, The Bridge River assemblage in the Meager Mountain volcanic complex, southwestern British Columbia, Current Research, Part E, Geol. Survey Can., Paper 90-1E: 153-157, 1990.  Stormer, J.C. , and J. Nicholls, XLFRAC: a program for the interactive testing of magmatic differentiation models, Comput. Geosci., 4, 143-159, 1978. Stout, M.Z., and J. Nicholls, Mineralogy and petrology of Quaternary lavas from the Snake River Plain, Idaho, Can. J. Earth Sci., 14, 2140-2156, 1977.  86  Wright, T.L., and P.C. Doherty, A linear programming and least square computer method for solving petrologic mixing problems, Geol. Soc. America. Bull., v. 81, 1995-2008, 1970.  87  APPENDICES  Al Linear Algebra Theorem of SVD  SVD is based on the following theorem of linear algebra: Any MxN matrix (MxN) A can be written as the product of an MxN column-orthogonal matrix U, an NxN diagonal matrix W' with positive or zero element, and the transpose of an NxN orthogonal matrix V ( Press et al., 1987). This theorem can simply be expressed by:  A =- U•W•VT. -  88  A2 Detailed Petrographic Descriptions The petrography of each MM sample is studied and the detail information is described below. A2.1 Pliocene assemblage (PLO)  PLO rock members are: MM-24 and MM36.  The Pliocene Assemblage was formed during the late Tertiary volcanic episode. As the oldest unit in this volcanic area, contacting with the basement metamorphic unit directly, it is almost covered by the younger units. Only four outcrop traces of flows, intrusions and breccias are exposed, which is mostly appears at the southern and southwestern edge of the map area, except one plagioclase andesite flow is exposed in the western Affliction Creek. At least three flow are recognized within this assemblage. The stratigraphy can be expressed as following. The bottom layer of medium grey-green aphanitic flows and breccias is overlaid by the volcanic breccia with dominant pluton and some aphanitic volcanic clasts which is as thick as up to 300m. The top two layers are plagioclase andesite flows and breccias, plagioclase, quartz and hornblende dacite flows and breccias in decreasing age order (Read, 1977). Samples yield the age about 1.9 to 2.2 Ma with an error of 0.2Ma by K-Ar dating (Green, Armstrong et al., 1988). Two samples were collected from the top two layers of this unit. One is from the plagioclase andesite flows (2.0+-0.1Ma), another from the quartz, plagioclase and hornblende dacite flow (1.9 +- 0.2Ma). Although only two samples were collected from this area, they are quite different from the mineralogy and petrology. MM-36 came from the andesite unit. The detailed information about the thin sections can be inferred from the appendix A4.1. The most distinguish feature  89  is: MM-24 has the quartz and biotite phenocrysts and MM-36 doesn't. They are common in that both have the plagioclase and hornblende phenocrysts, clinopyroxene microphenocryst. Sample number: MM-24 Thin section: Plagioclase phenocryst, quartz xenocryst, clinopyroxene, hornblende and biotite phenocrysts, minor apatite and zircon microphenocrysts, quartz glomerocryst are carried in the groundmass of plagioclase, hornblende, opaques and glass. Anhedral to subhedral plagioclase is zoned with the size from 2mm to 0.5mm. Most of the pl grains have 0.01mm wide remelted fringe and very thin overgrowth rim, except one or two small pl phenocrysts without the overgrowth rim. All the quartz grains are embayed with the fractured surface, and the melt corroded the quartz from the cracks. The disequilibrium relation between the quartz and the melt and the absence of the quartz in the groundmass infer that the quartz grains are the xenocrysts. The biotites were broken down to plagioclases and Fe oxides. The hornblende both as the twinned phenocrystals and in the groundmass were altered with the opaque rim. Some of the hornblende microphenocrystals were totally replaced by the opaques. The number of the clinopyroxene is a small amount, only two were found in one thin section, but the size is about 0.15mm, not so small. There is not any reaction rim and alteration phenomenon. One quartz glomerocryst about 0.07mm in diameter with minor hornblende is carried in the groundmass, the hornblende was almost altered to the opaques. One zircon microphencryst was noticed on the groundmass with the unusual interference color. Mineral classification: Mineral^Mode^Name Plagioclase^5%^Andesite Quartz^2% Hornblende^1% Biotite^1% 90  Clinopyroxene^1% Groundmass^90% Sample numbler: MM-36 Thin section: Plagioclase, hornblende, clinopyroxene and opaques phenocrysts, plagioclase, clinopyroxene and opaques glomerocryst are carried in the groundmass of plagioclase, hornblende, orthopyroxene, clinopyroxene, opaques and glass. Plagioclase varies from 0.15mm to 0.03mm in size. Two types of plagioclase exist in this thin section. One type with the remelted fringe and the overgrowth rim, while the other doesn't have, but both types are strongly zoned. The hornblende appears the light yellow to deep reddish brown pleochroism. All the euhedral, subheral and anhedral hornblende grains are fresh without any reaction rim and alteration surface. While the orthopyroxene has the hornblende coronas. The irregular habit of the opaques and their coexisting with the ferromagnesian mineral infer that they came from the alteration. The clot comprises the orthopyroxene, opaque and plagioclase. The orthopyroxene has the hornblende reaction rim, the opaques keep the pyroxene shape wirh the reddish brown hornblende rim. These facts suggests that the clots are produced by orthopyroxene breakdown. Mineral classification: Mineral^Mode^Name Plagioclase^10%^Dacite Orthopyroxene^2% Hornblende^3% Opaques^1% Clot^0.5% Groundmass^87.5%  91  A2.2 Pylon Assemblage (PYA) PYA rock members are: MM-20, MM-21, MM-25, MM-28 and MM-35.  Sample number: MM-20 Thin section: Plagioclase, hornblende, orthopyroxene, biotite and opaques phenocrysts, plagioclase and hornblende glomerocryst are carried in the groundmass of plagioclase, hornblende and glass. Subhedral to anhedral plagioclase varies from 0.2mm to 0.05mm in size, slightly zoned with the interference color high to first-order yellow. Glass inclusion is rich in the plgioclase grains. Euhedral to subhedral hornblende has the light greenish yellow to deep brown pleochroism with the apatite inclusion. Minor anhedral biotite was slightly remelted. Microphenocrysts of orthopyroxene were found in the groundmass. Most of the opaques have the irregular shape, a few grains occur as the nice square, which are possible magnetite. The groundmass is glassy and rich in plagioclase. The plagioclase and hornblende in the glomerocryst are the same as the phenocrysts. And there isn't any reaction or alteration relation is noticed. Therefore this crystal clot is formed by the minerals growing together. Mineral classification: Mineral^Mode^Name Plagioclase^10%^Andesite Hornblende^3% orthopyroxene^1% opaque^1% glomerocryst^0.5% biotite^0.5% Groundmass^84%  92  Sample number: MM-21 Thin section: Plagioclase, hornblende, orthopyroxene and opaque phenocrysts, plagioclase, orthopyroxene, hornblende, biotite and opaque glomerocrysts are carried in the groundmass of plagioclase, orthopyroxene, opaque and glass. The subhedral to anhedral plagioclases are common zoned, with the glass liquid inclusion. The rim is round because of the remelted. The hornblende has the light yellow green to deep green pleochroism. Most grains are andedral with the remelted rim. A few have the orthopyroxene coronas. All the orthopyroxene are anhedral and small in size about 0.02mm. Two glomerocrysts are noticed, one with biotite. Both are common with hornblende, plagioclase, opaque and orhtopyroxene. The hronblende is the same in clots and as phenocrysts. Therefore these two clots are the production of the hornblende breakdown. Mineral classification: Mineral^Mode^Name Plagioclase^10%^A ndesite Hornblende^3% Orthopyroxene^1% Opaque^0.5% Glomerocryst^1.5% Groundmass^84% Sample number: MM-25 Thin section: Plagioclase, hornblende, orthopyroxene, clinopyroxene and opaques phenocrysts, plagioclase, orthopyroxene, clinopyroxene and opaques glomerocryst are carried in a glassy groundmass of plagioclase and opaques. Subhedral to anhedral plagioclase zonedly grew with the needle like apatite, second-order interferene color clnopyroxene, first-order interference orthopyroxene, and higher-order white interference calcite. The melt embayed the plagioclase from the (010) 93  direction. The hornblende does't appear as plgioclase inclusion, but euhedral to subhedral phenocryst with size about 0.05mm. With the light yellow to deep brown pleochroism in common, all the hornblende grains were oxided with the Fe-oxide rim. Some of the small grains were totally replaced by the opaques. The euhedral twinned clinopyroxene and subhedral elongate clinopyroxene almost have the equal amount as the orthopyroxene. The minerals in the glomerocryst are the same as phenocrysts, but disequilibrium is obvious existing within this crystal clot. This infer that the glomerocryst is probably the production of settlement of plagioclase and pyroxene at the bottom of the magma chamber. Mineral classification: Mineral^Mode^Name Plagioclase^15 %^Andesite Clinopyroxene^3% Orthopyroxene^3% Hornblende^2% Opaques^1% Groundmass^76% Sample number: MM-28 Thin section: Plagioclase, hornblende, orthopyroxene and opaques phenocrysts, clinopyroxene and apatite microphenocrysts, plagioclase and orthopyroxene glomerocryst, plagioclase and opaques glomerocryst are carried in the groundmass of plagioclase, orthopyroxene, clinopyroxene, opaques and glass. The plagioclase has the minor apatite and pyroxene inclusion, showing sieved structure. The sieved plagioclase varies with the remelted fringe wide from 0.1mm to 0.01mm. Hornblende is altered to opaques. Euhedral to subhedral orthopyroxene is melted to plagioclase. Mineral classification: 94  ^ Name Mineral^Mode ^ Andesite Plagioclase^10% Hornblende^1% Orthopyroxene^2% Opaques^1% Groundmass^86% Sample number: MM-35 Thin section: Plagioclase, hornblende and orthoroxene phenocrysts, plagioclase, hornblende and orthopyroxene glomerocryst are carried in the groundmass of plagioclase, hornblende, orthopyroxene and glass. This sample is very vesicular. Plagioclase grew in the centre of a vesicle. Hornblende has the light yellow to deep red pleochroism. All the ferromagnesian mineral, like hornblende and orthopyroxene were altered to opaques. Mineral classification: Mineral^Mode^Name Plagioclase^15 %^Andesite Hornblende^3% Orthopyroxene^2% Opaques^1% Groundmass^79% A2.3 Capricorn Assemblage (CPA)  CPA rock members are: MM-34 and MM-38. Sample number: MM-34 95  Thin section: Plagioclase, quartz, biotite and hornblende pheonocrysts, apatite microphenocryst are carried in the groundmass of plagioclase, hornblende, orthopyroxene, opaques and glass. The plagioclase shows two different groups. One group has very thin multiple twin lathe. Another group does't have the obvious mutiple twin, but nice zoned. Two types of plagioclase have the sieved texture, shown by the remelted fringe and overgrowth rim. Nice glass inclusion is found in a few plagioclase grains. One plagioclase ia found including zircon microphenocryst. The quartz is embayed and corroded by the fluid and magmatic liquid. Fratured surface hint the strong strain. Small quartz piece is also found in the groundmass. This infer that quartz is one of the early crystallizing mineral phase. The hornblende and biotite both are altered to opaques. The microphenocryst and groundmass orthopyroxene is also oxidized at the rim. Mineral classification: Mineral^Mode^Name Plagioclase^7%^Rhyodacite Quartz^2% Hornblende^1% Biotite^1% groundmass^89% Sample number: MM-38 Thin section: Plagioclase and hornblende phenocrysts are carried in the highly vesicular groundmass of plagioclase. Subhedral plagioclase has an average size of 0.03mm with the largest up to 0.1mm. This sample is hightly vesicular. Plagioclase is slightly zoned, embayed at the rim. Sieved texture is obviously in this sample, early disequilibrium shown as the irregular remelted inside zone, later stable condition brought an overgrowth rim. The yellow green hornblende was resorbed, with an opaque rim. The lackage of hornblende in the groundmass suggests that 96  the hornblende is possibly xenocryst. The crystallizing sequence is probably as plagioclase, then amplibole followed by magnetite. Mineral classification: Mineral^Mode^Name Plagioclase^3%^Rhyolite Hornblende^0.5% Vesicle^10% Groundmass^86% A2.4 Devastation Assemblage (DVA)  DVA rock members are: MM-27 and MM-33.  Sample number: MM-27 Thin Section: Plagioclase, hornblende and orthopyroxene phenocrysts, hornblende, plagioclase and opaques glomerocryst are carried in the groundmass of plagioclase, orthopyroxene, clinopyroxene and glass. The subhedral plagioclase has the size from 0.2mm to 005mm, slightly zoned, sieved, embayed. The hornblende is reacting to plagioclase and opaque with liquid, While it keeps the hornblende's shape. Orthopyroxene has slight pleochroism, from pale to pale pink, very low interference color. Mineral classification: Mineral^Mode^Name Plagioclase^15% Hornblende^3% Orthopyroxene^2% Groundmass^80% 97  Sample number: MM-33 Thin section: Plgioclase, quartz and biotite phenocrysts are carried in the groundmass of plagioclase, biotite and glass. Two groups of plagioclase appear. One show the sieved texture, with the remelted fringe and overgrowth rim, slightly zoned, slightly twinned. Another shows the thin multiple twinning lathe, and zoned, but without the remelted fringe and overgrowth rim. Biotite is being broken down by the melt. Subhedral to anhedral quartz is remelted, embayed, by the disequilibrium liquid. Mineral classification: Mineral^Mode^Name Plagioclase^12% Quartz^2% Biotite^1% Groundmass^85 % A2.5 Plinth Assemblage (PLA)  PLA rock members are: MM-26, MM-29, MM-30, MM-39 and MM-40. Sample number: MM-26 Thin section: Plagioclase, quartz, hornblende and biotite phenocrysts, apatite microphenocryst are carried in the groundmass of plagioclase, hornblende, clinopyroxene, opaques and glass. Plagioclase is zoned and the grains of plagioclase tend to attach together by the sysnesis process. The fractured quartz is embayed, corroded and even with the clinopyroxene coronas. Apatite inclusions are visible in some grains of quartz. Anhedral biotite has the reaction rim of 98  hornblende and minor plagioclase. Hornblende is replacing the biotite. Microphenocrysts of clinopyroxene and apatite are rich in the groundmass. Mineral classification: Mineral^Mode^Name Plagioclase^15%^Rhyodacite Quartz^5% Hornblende^1% Biotite^3% Clinopyroxene^0.5% Apatite^0.1 % Groundmass^75.4% Sample number: MM-29 Thin section: Plagioclase, quartz, biotite and hornblende phenocrysts, hornblende, plagioclase and clinopyroxene microphenocrysts are carried in the groundmass of plagioclase, hornblende and glass. There are two types of plagioclase. One is zoned with the thin lathe of multiple twinning. Another is zoned but with the coarse core. The quartz has the fractured surface, embayed and corroded by the melt. Euhedral to subhedral biotite has the size about 10.5m, some grains have the hornblende reaction rim. Hornblende has a much smaller size relatively to biotite. Twinned clinopyroxene microphenocryst scattered in the groundmass. Mineral classification: Mineral^Mode^Name Plagioclase^20%^Rhyodacite Quartz^3% Biotite^5% Hornblende^2% 99  Clinopyroxene^1% Groundmass^69% Sample number: MM-30 Thin section: Plagioclase, quartz, biotite and hornblende phenocrysts, clinopyroxene and orthopyroxene microphenocrysts are carried in the groundmass of plagioclase, hornblende and glass. The plagioclase has the sieved texture, shown as the remelted fringe and overgrowth rim. The quartz has the clinopyroxene reaction rim. Hornblende is replacing the biotite, shown as the biotite has the hornblende reaction rim. Mineral classification: Mineral^Mode^Name Plagioclase^12%^Andesite Quartz^7% Biotite^5% Hornblende^2.5% Groundmass^73.5% Sample number: MM-39 Thin section: Plagioclase, quartz and biotite phenocrysts, clinopyroxene, opaque and zircon microphenocrysts are carried in the groundmass of plagioclase, opaques and glass. Apatite inclusion in plagioclase is visible. The groundmass is in red representing the oxidized status of this sample. Biotite is breaking down to opaques. Mineral classification: Mineral^Mode^Name Plagioclase^12%^Rhyodacite 100  Quartz^2% Biotite^3% Clinopyroxene^1% Opaques^1% Groundmass^81 % Sample number: MM-40 Thin section: Plagioclase, quartz and biotite phenocrysts, hornblende microphenocryst are carried in the groundmass of plagioclase, hornblende and glass. There are very small amount of pyroxene, mostly clinopyroxene on the groundmass. Biotite is breaking down to opaque. Both plagioclase and quartz are remelted and with the overgrowth rim. The quartz is regarded as xenocrystal. Mineral classification: Mineral^Mode^Name Plagioclase^17%^Rhyodacite Biotite^4% Quartz^4% Hornblende^1% Groundmass^74% A2.6 The Mosaic Assemblage (MOA)  MOA rock members are: MM-17, MM-18, MM-19 and MM-37. Sample number: MM-17 Thin section: (located at north exposure) 101  Plagioclase and olivine phenocrysts are carried in the groundmass of plagioclase, olivine and titanaugite. Plagioclase is slightly zoned. Titanaugite fill between the plagioclase and olivine in the intergranular texture. Plagioclase is slightly zoned. Mineral classification: Mineral^Mode^Name Plagioclase^15 %^Alkali olivine basalt Olivine^15 % Titanaugite^10% Groundmass^60% Sample number: MM-18 Thin section: (located at south exposure) Big olivine and plagioclase phenocrysts are carried in the groundmass of plagioclase, olivine and clinopyroxene. Olivine is slightly zoned and with alteration rim. Mineral classification: Mineral^Mode^Name Plagioclase^15 %^Alkali olivine basalt Olivine^15% Groundmass^70% Sample number: MM-19 Thin section: (located at south exposure) Plagioclase and olivine phenocrysts are carried in the groundmass of plagioclase, clinopyroxene and olivine. Olivine has the clean surface and clear rim. Magnetite appears on the surface of olivine. Mineral classification: Mineral^Mode^Name 102  Olivine^15%  ^  Alkali olivine basalt  Plagioclase^10% Ground mass^75% Sample number: MM-37 Thin section: (located at west exposure) Plagioclase, olivine and clinopyroxene phenocrysts are carried in the groundmass of plagioclase, olivine and clinopyroxene. Minor anorthoclase or sanidine with melting rim occur in the groundmass. Plagioclase has the melting overgrowth surface. Mineral classification: Mineral^Mode^Name Plagioclase^10%^Hawaiite Olivine^5% Augite^10% groundmass^75 %  103  A3 Whole-Rock Chemical Analysis The operation conditions for the XRF analysis of major, minor and trace elements are in Table A3.1 and A3.2. Table 3.3 is the CIPW calculations for the four rocks in Mosaic assemblage in order to present the influence of iron status on CIPW results.  104  A4 Electron Microprobe Analysis Mineral chemistry was studied on the scanning electrical microprobe. The chemical components of phenocrystals provide the essential data for mass balance calculation. The compositions of plagioclase, biotite, amphibole and pyroxene were measured and the analytical results are presented in Table 4.1, 4.2, 4.3 and 4.4.  105  Table A3-1. XRF operating conditions for the analysis of major and minor elements.  Fe pk Fe bkg Ti pk Ti bkg Ca pk Ca bkg K pk K bkg Si pk Si bkg Line^Ica^--^ka^--^ka^--^ka^--^ka^-20^57.49 56.00^86.12 91.00^113.12 110.40^136.70 132.15 109.24 113.3 0 Target^Cr^Cr^Cr^Cr^Cr Cr^Cr^Cr^Cr^Cr Crystal^LIF200 LIF200 LIF200 LIF200 LIF20 LIF200 LIF200 LIF200 PET PET 0 Icv/mA^60/40 60/40^60/40 60/40 60/32 60/32 ^60/40 60/40^60/40 60/40 Collimator^C^C^C^C^C^C^C^C^C^C Counter^F^F^F^F^F^F^F^F^F^F Vacuum^On On^On On^On On^On^On^On On Gain^128^128^128^128^128^128^128^128^128^128 Counterkv^830x2 830x2^830x2 830x2 830x2 830x2^830x2 830x2 830x2 830x2 Lower level^200^200^200^200^200 200^200^200^200 200 Window^800^800^800^800^800 800^800^800^800 800 Count time^lOsec lOsec^lOsec lOsec^lOsec lOsec^lOsec lOsec^100sec lOsec  Element^Al pk Al bkg P pk P bkg Mg pk Mg bkg Na pk Na bkg Mn pk Mnbkg Line^ka^--^ka^--^ka^--^ka^--^ka^-20^145.02 139.00^89.53 92.60 45.15 44.00^55.20 53.30^62.97 65.00 Target^Cr^Cr^Cr^Cr^Cr Cr^Cr^Cr^Mo Mo Crystal^PET PET^PET PET TIAP TIAP^TIAP TIAP LIF20 LIF200 0 kv/mA^60/40 60/40^60/45 60/45 60/45 60/45^60/45 60/45^60/45 60/45 Collimator^C^C^C^C^C^C^C^C^C^C Counter^F^F^F^F^F^F^F^F^F^F Vacuum^On^On^On On^On On^On^On^On On Gain^128^128^128^128^128^128^128^128^128^128 Counter kv^830x2 830x2^830x2 830x2 830x2 830x2^830x2 830x2 830x2 830x2 Lower level^200^200^200^200^200 200^200^200^200 200 Window^800^800^800 800^800 800^800^800^800 800 Count time^100sec lOsec^100sec lOsec^100sec lOsec^100sec lOsec^100sec lOsec  * The XRF analyses were done on XRF spectrometer PW 1410 in the Department of Geological Sciences of UBC.  106  Table A3-2 XRF operating conditions for the analysis of trace elements.  Element  Nb pk  Nb -bkg  Nb+bkg  Y pk  Y -bkg^Y +bkg  Line 20 Target Crystal kv/mA Collimator Counter Vacuum Gain Counter kv Lower level Window Count time  La 21.40 Cr LIF200 60/40 F S.I.  21.10 Cr LIF200 60/40 F S.I.  21.70 Cr LIF200 60/40 F S.I.  On  La 23.80 Cr LIF200 60/40 F S.I.  On  -23.50^24.10 Cr^Cr LIF200^LIF200 60/40^60/40 F^F S.I.^S.I.  On  128 1044 280 420 40sec  128 1044 280 420 40sec  On  On^On  128 1044 280 420 40sec  128 1044 280 420 40sec  128^128 1044^1044 280^280 420^420 40sec^40sec  Element  Zr pk  Sr pk  Rb pk  Zr,Sr,Rb -bkg Zr,Sr,Rb+bkg  Line 20 Target Crystal kv/mA Collimator Counter Vacuum Gain Counter kv Lower level Window Count time  La 22.51 Cr LIF200 60/40 F S.I.  La 25.15 Cr LIF200 60/40 F S.I.  La 26.60 Cr LIF200 60/40 F S.I.  19.50 Cr LIF200 60/40 F S.I.  On  On  24.55 Cr LIF200 60/40 F S.I.  On  On  128 1044 280 420 40sec  128 1044 280 420 40sec  On  128 1044 280 420 40sec  128 1044 280 420 40sec  128 1044 280 420 40sec  * The XRF analyses were carried on the XRF spectrometer in the Department of Geological Sciences of UBC.  Table A3-3 Effect of Fe 2 +/Fe3 + ratio on the computed CIPW norm for the rocks of Mosaic assemblage  Fe state Sample Ne Or P1 Di Hy 01 II Hm Tn Ap Mt  17  Fe total as Fe 7 0 / 18 19  37  17  Fe total as FeO 18 19  37  0.00 0.00 0.00 0.00 1.90 1.21 1.48 4.83 4.32 4.57 4.55 7.76 4.32 4.57 4.55 7.76 66.55 57.57 56.52 60.83 63.39 55.54 54.06 52.77 7.56 11.66 10.80 12.71 10.79 15.14 14.46 16.46 8.11 12.11 11.97 0.11 0.00 0.00 0.00 0.00 3.49 2.99 4.63 8.00 17.03 20.97 22.79 15.03 0.23 0.25 0.25 0.22 1.85 1.99 2.08 2.09 6.59 7.67 7.95 6.50 0.00 0.00 0.00 0.00 2.42 2.61 2.74 2.81 0.00 0.00 0.00 0.00 0.72 0.56 0.59 1.21 0.72 0.56 0.59 1.06 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00  108  Fe total as Fe.,0 / + FeO 17 18 19 37 0.00 0.00 0.00 1.87 4.32 4.57 4.55 7.76 66.55 57.57 56.52 57.71 10.79 15.14 14.46 16.46 0.04 2.02 1.08 0.00 13.18 15.52 18.22 9.10 1.85 1.99 2.08 2.09 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.72 0.56 0.56 1.06 2.63 2.55 2.51 3.95  Table A4-1 Electron microprobe analyses of PLA plagioclases.  Sample Oxides  MM-26A core  MM-26A rim  MM-26B core  MM-26B core  MM-26B rim  MM-26 average  SiO 2 Al 2 03 FeO MgO Na2 O K2 O CaO  61.05 24.40 0.15 0.01 7.55 0.51 6.28  61.17 24.31 0.27 0.01 7.48 0.59 5.98  63.68 22.78 0.20 0.00 8.32 0.85 4.13  63.57 22.86 0.25 0.00 8.14 0.79 4.40  63.36 22.94 0.22 0.00 8.11 0.79 4.69  62.57 23.46 0.22 0.00 7.92 0.71 5.10  Total  99.95  99.81  99.96  100.00  100.11  99.97  Si 4 + A1 3 +  2.7174 1.2798 3.9972  2.8028 1.1956 3.9987  2.7750 1.2262 4.0012  K4 Na+ Ca2 + mg2+  0.0291 0.6514 0.2995 0.0008 0.0054 0.9862  0.0337 0.6463 0.2856 0.0004 0.0099 0.9759  0.0478 0.7132 0.1958 0.0000 0.0074 0.9642  0.0444 0.6976 0.2085 0.0003 0.0091 0.9600  0.0445 0.6958 0.2225 0.0000 0.0082 0.9710  0.0399 0.6810 0.2423 0.0003 0.0080 0.9715  K Na Ca  0.0291 0.6541 0.2995 0.9801  0.0337 0.6463 0.2856 0.9656  End Members 0.0478 0.7132 0.1958 0.9568  0.0444 0.6976 0.2085 0.9506  0.0445 0.6958 0.2225 0.9628  0.0399 0.6810 0.2423 0.9632  Others  0.0062  0.0104  0.0074  0.0094  0.0082  0.0083  Fe2 +  Structural Formula (Based on 8 Oxygens) 2.7249 2.8173 2.8119 1.2763 1.1879 1.1917 4.0012 4.0052 4.0036  109  Table A4.1 continued.  Sample Oxides  MM-29A 1  MM-29A 2  MM-29A 3  MM-29A 4  MM-29A 5  MM-29C core  MM-29C rim  MM-29B 1  MM-29B 2  MM-29B 3  MM-29B 4  MM-29B 5  MM-29 average  Si0 2 Al2 03 FeO MgO Na2 O K2 0 CaO  57.45 27.17 0.22 0.01 5.97 0.35 9.12  57.36 26.79 0.18 0.01 6.05 0.28 8.83  60.26 24.76 0.31 0.01 7.14 0.45 6.58  61.37 23.94 0.16 0.00 7.46 0.64 5.79  62.66 23.41 0.15 0.00 8.10 0.73 4.94  61.99 23.25 0.22 0.00 8.00 0.72 4.84  63.10 22.85 0.15 0.00 8.22 0.71 4.39  62.34 23.44 0.23 0.00 7.92 0.82 4.93  62.56 23.28 0.14 0.01 8.05 0.71 4.87  61.34 23.54 0.09 0.01 7.70 0.74 4.97  63.43 22.39 0.17 0.01 8.47 0.95 3.76  63.49 22.26 0.17 0.00 8.53 0.94 3.67  61.44 23.92 0.18 0.01 7.63 0.67 5.56  Total  100.29  99.49  99.51  99.38  99.99  99.02  99.43  99.68  99.62  98.38  99.19  99.06  99.42  Si 4 + A1 3 +  2.5690 1.4319 4.0009  2.5818 1.4211 4.0029  2.6957 1.3055 4.0012  2.7420 1.2608 4.0027  Structural Formula (Based on 8 Oxygens) 2.7780 2.7756 2.8070 2.7739 1.2235 1.2269 1.1978 1.2291 4.0015 4.0025 4.0048 4.0030  2.7822 1.2205 4.0027  2.7625 1.2394 4.0119  2.8279 1.1764 4.0043  2.8332 1.1708 4.0039  K+ Na+ Ca2 + mg2+  0.0198 0.5179 0.4369 0.0004 0.0084 0.9834  0.0164 0.5279 0.4256 0.0008 0.0066 0.9770  0.0259 0.6196 0.3156 0.0008 0.0116 0.9734  0.0368 0.6467 0.2773 0.0000 0.0061 0.9669  0.0414 0.6963 0.2346 0.0002 0.0054 0.9780  0.0409 0.6946 0.2323 0.0003 0.0083 0.9765  0.0402 0.7090 0.2091 0.0000 0.0059 0.9642  0.0465 0.6835 0.2350 0.0000 0.0087 0.9737  0.0404 0.6938 0.2320 0.0008 0.0053 0.9723  0.0428 0.6723 0.2397 0.0004 0.0034 0.9587  0.0542 0.7319 0.1798 0.0009 0.0062 0.9730  0.0538 0.7382 0.1754 0.0000 0.0062 0.9737  Na Ca  0.0198 0.5179 0.4369 0.9746  0.0164 0.5276 0.4256 0.9696  0.0259 0.6196 0.3156 0.9611  0.0368 0.6467 0.2773 0.9608  0.0414 0.6963 0.2346 0.9723  End Members 0.0409 0.0402 0.6946 0.7090 0.2323 0.2091 0.9583 0.9678  0.0465 0.6835 0.2350 0.9650  0.0404 0.6938 0.2320 0.9662  0.0428 0.6723 0.2397 0.9548  0.0542 0.7319 0.1798 0.9659  0.0538 0.7382 0.1754 0.9674  Others  0.0089  0.0074  0.0123  0.0061  0.0057  0.0086  0.0059  0.0087  0.0060  0.0039  0.0071  0.0062  Fe2 +  K  110  Table A4.1 continued.  Sample Oxides  MM-30B MM-30B MM-30B MM-30B MM-30B MM-30B MM-30B MM-30B MM-30B MM-30B MM-30B MM-30A MM-30AMM-30AMM-30C MM-30C MM-30 middle^rim^core^rim^average 1^2^3^4^5^6^7^8^9^10^11^core  SiO 2 Al 2 03 FeO MgO Na2 O K2 O CaO  63.77 22.64 0.24 0.01 8.18 0.82 4.42  54.88 28.39 0.18 0.00 5.03 0.28 10.99  56.42 27.62 0.18 0.00 5.62 0.251 9.68  56.48 27.43 0.20 0.02 5.74 0.29 9.76  60.86 24.71 0.22 0.00 7.34 0.53 6.11  62.87 23.21 0.19 0.02 7.92 0.73 4.96  56.83 26.96 0.15 0.01 5.97 0.30 8.63  62.48 23.51 0.22 0.01 7.77 0.75 5.14  59.35 25.73 0.22 0.02 6.72 0.45 7.77  61.86 24.07 0.12 0.00 7.48 0.58 5.60  63.74 22.82 0.22 0.02 8.15 0.96 4.06  62.67 23.14 0.17 0.01 8.17 0.82 4.74  62.67 23.43 0.23 0.02 8.15 0.76 4.56  59.52 24.99 0.45 0.05 5.61 0.51 8.35  62.97 22.93 0.13 0.00 7.98 0.89 4.62  62.60 23.48 0.17 0.00 7.74 0.69 4.92  60.62 24.69 0.21 0.01 7.10 0.60 6.52  Total  100.08  99.69  99.75  99.92  99.78  99.91  98.85  99.88  100.24  99.71  99.96  99.71  99.83  99.48  99.52  99.60  99.74  Si 4 + A1 3 +  2.8192 1.1798 3.9991  2.4822 1.5135 3.9957  2.5389 1.4646 4.0035  2.5405 1.4541 3.9946  2.7112 1.2973 4.0085  2.7878 1.2128 4.0006  2.5732 1.4386 4.0118  2.7733 1.2300 4.0032  2.6450 1.3515 3.9965  2.7494 1.2610 4.0104  2.8191 1.1893 4.0083  2.7868 1.2125 3.9994  2.7814 1.2255 4.0070  2.6696 1.3209 3.9904  2.8013 1.2023 4.0036  2.7812 1.2295 4.0107  2.7045 1.2982 4.0027  K+ Na+ Ca2 + Mg2 + Fe2 +  0.0462 0.7013 0.2093 0.0003 0.0089 0.9661  0.0159 0.4412 0.5301 0.0002 0.0066 0.9941  0.0144 0.4901 0.4665 0.0000 0.0066 0.9775  0.0167 0.5006 0.4703 0.0016 0.0077 0.9968  0.0302 0.6342 0.2915 0.0000 0.0081 0.9641  0.0413 0.6807 0.2358 0.0015 0.0071 0.9664  Structural Formula (Based on 8 Oxygens) 0.0176 0.0422 0.0254 0.0329 0.0544 0.5241 0.6691 0.5802 0.6445 0.6985 0.4186 0.2445 0.3710 0.2667 0.1926 0.0007 0.0003 0.0010 0.0000 0.0010 0.0058 0.0084 0.0081 0.0046 0.0080 0.9668 0.9645 0.9857 0.9486 0.9546  0.0465 0.7042 0.2258 0.0004 0.0062 0.9831  0.0428 0.7014 0.2169 0.0016 0.0085 0.9712  0.0292 0.4881 0.4013 0.0032 0.0167 0.9385  0.0503 0.6881 0.2203 0.0000 0.0048 0.9635  0.0393 0.6669 0.2344 0.0000 0.0062 0.9467  0.0342 0.6139 0.3115 0.0007 0.0077 0.9680  K Na Ca  0.0462 0.7013 0.2093 0.9568  0.0159 0.4412 0.5301 0.9873  0.1444 0.4901 0.4663 0.9710  0.0167 0.5006 0.4703 0.9876  0.0302 0.6342 0.2915 0.9560  0.0413 0.6807 0.2358 0.9578  0.0176 0.5241 0.4186 0.9603  End Members 0.0422 0.0254 0.6691 0.5802 0.2445 0.3710 0.9558 0.9766  0.7042 0.2258 0.9765  0.7014 0.2169 0.9611  0.4881 0.4013 0.9186  0.6881 0.2203 0.9586  0.6669 0.2344 0.9405  0.6139 0.3115 0.9596  0.0067  0.0101  0.0200  0.0048  0.0062  0.0084  Others  0.0093  0.0068  0.0066  0.0093  0.0081  0.0086  0.0065  0.0087  111  0.0091  0.0329 0.6445 0.2667 0.9440  0.0544 0.6985 0.1926 0.9455  0.0046  0.0090  Table A4.1 continued.  Sample Oxides  MM-39B 1  MM-39B 2  MM-39B 3  MM-39B 4  MM-39B 5  MM-39B 6  MM-39B 7  MM-39B 8  MM-39C 1  MM-39C 2  MM-39C 3  MM-39C MM-39C 4^5  Si0 2 Al 2 0 3 FeO MgO Na2 O K2 0 CaO  62.37 23.51 0.20 0.00 7.98 0.65 5.18  62.08 24.16 0.12 0.01 7.84 0.67 5.43  61.30 24.04 0.17 0.00 7.65 0.76 5.10  61.35 24.73 0.267 0.00 7.89 0.84 4.62  61.03 24.32 0.19 0.00 7.81 0.81 4.92  60.54 24.48 0.28 0.00 7.52 0.72 5.44  59.46 25.47 0.23 0.00 7.55 0.80 5.32  62.41 23.64 0.15 0.00 7.64 1.01 4.80  62.91 23.33 0.19 0.00 7.94 0.79 4.74  62.43 23.27 0.14 0.00 8.05 0.79 4.67  62.80 23.27 0.15 0.00 8.25 0.77 4.51  63.23 22.84 0.18 0.00 8.40 0.82 4.08  61.63 24.18 0.15 0.00 8.08 0.81 4.31  Total  99.90  100.32  99.01  99.70  99.08  98.98  98.84  99.66  99.89  99.34  99.76  99.56  99.17  Structural Formula (Based on 8 Oxygens) 2.7343 2.7185 2.6775 2.7747 1.2388 1.2840 1.2955 1.3515 4.0183 4.0140 4.0290 4.0135  2.7885 1.2187 4.0071  2.7837 1.2228 4.0065  2.7882 1.2179 4.0060  2.8097 1.1963 4.0060  2.7529 1.2730 4.0259  Si 4 + A1 3 +  2.7693 1.2304 3.9996  2.7465 1.2596 4.0061  2.7458 1.2681 4.0149  2.7301 1.2968 4.0269  K+ Na+ Ca2 + m g2 + pe2+  0.0366 0.6874 0.2465 0.0000 0.0076 0.9781  0.0379 0.6727 0.2574 0.0004 0.0046 0.9730  0.0433 0.6644 0.2448 0.0000 0.0063 0.9587  0.0479 0.6811 0.2204 0.0000 0.0100 0.9593  0.0465 0.6783 0.2361 0.0000 0.0070 0.9680  0.0410 0.6547 0.2619 0.0000 0.0103 0.9679  0.0462 0.6591 0.2565 0.0003 0.0086 0.9708  0.0576 0.6585 0.2289 0.0001 0.0057 0.9505  0.0448 0.6825 0.2249 0.0000 0.0070 0.9591  0.0452 0.6960 0.2230 0.0000 0.0053 0.9694  0.0436 0.7101 0.2147 0.0000 0.0057 0.9742  0.0466 0.7241 0.1943 0.0000 0.0067 0.9717  0.0461 0.6998 0.2064 0.0003 0.0055 0.9580  K Na Ca  0.0366 0.6874 0.2465 0.9705  0.0379 0.6727 0.2574 0.9679  0.0433 0.6644 0.2448 0.9525  0.0479 0.6811 0.2204 0.9493  0.0465 0.6783 0.2361 0.9610  End Members 0.0410 0.0462 0.6591 0.6547 0.2619 0.2565 0.9575 0.9618  0.0573 0.6585 0.2289 0.9447  0.0448 0.6825 0.2249 0.9521  0.0452 0.6960 0.2230 0.9641  0.0436 0.7101 0.2147 0.9684  0.0466 0.7241 0.1943 0.9650  0.0461 0.6998 0.2064 0.9522  Others  0.0076  0.0051  0.0063  0.0100  0.0070  0.0103  0.0059  0.0070  0.0053  0.0057  0.0067  0.0058  112  0.0090  Table A4.1 continued.  Sample Oxides  MM-39C 6  MM-39C MM-39C MM-39D MM-39D MM-39A MM-39A 7 8^core^rim^rim^rim  SiO 2 Al 2 0 3 FeO MgO Na2 O K2 O CaO  61.33 23.73 0.16 0.00 7.99 0.69 5.29  62.31 22.87 0.24 0.00 8.34 0.81 4.38  61.94 23.30 0.20 0.00 7.95 1.12 4.580  62.74 23.730 0.22 0.00 7.81 0.68 5.19  62.44 23.06 0.18 0.00 7.81 1.20 4.88  62.48^62.30 23.32^23.82 0.11^0.15 0.01^0.00 7.96^7.89 0.86^0.86 4.97^5.42  62.77 23.15 0.24 0.02 8.01 1.00 4.85  60.54 25.01 0.12 0.00 7.09 0.55 6.85  62.68 23.54 0.23 0.00 8.15 0.78 4.61  61.77 24.23 0.18 0.00 7.66 0.64 5.51  61.69 24.08 0.20 0.00 7.76 0.80 5.06  61.94 23.80 0.19 0.00 0.81 0.81 4.99  Total  99.14  98.93  99.11  100.37  99.56  99.72^100.45  100.04  100.16  100.00  100.00  99.59  99.61  Si 4 + Al 3 +  2.7480 1.2530 4.0010  2.7922 1.2079 4.0001  2.7750 1.2301 4.0051  2.7704 1.2348 4.0052  Structural Formula (Based on 8 Oxygens) 2.7852 2.7789^2.7559 2.7854 1.2122 1.2221^1.2421 1.2106 3.9974 4.0010^3.9980 3.9960  2.6913 1.3106 4.0019  2.7779 1.2297 4.0076  2.7413 1.2671 4.0084  2.7485 1.2643 4.0127  K+ Na+ Ca2 + mg2+  0.0397 0.6901 0.2540 0.0000 0.0059 0.9897  0.0462 0.7232 0.2103 0.0000 0.0089 0.9887  0.0644 0.6908 0.2198 0.0000 0.0077 0.9826  0.0383 0.6688 0.2457 0.0000 0.0080 0.9608  Structural Formula (Based on 8 Oxygens) 0.0680 0.0491^0.0485 0.0567 0.6753 0.6862^0.6770 0.6889 0.2333 0.2367^0.2570 0.2306 0.0000 0.0008^0.0000 0.0012 0.0065 0.0042^0.0056 0.0088 0.9831 0.9770^0.9880 0.9862  0.0312 0.6112 0.3263 0.0000 0.0044 0.9731  0.0444 0.7002 0.2190 0.0000 0.0087 0.9723  0.0364 0.6593 0.2620 0.0003 0.0068 0.9649  0.0457 0.6706 0.2414 0.0000 0.0075 0.9651  K Na Ca  0.0397 0.6901 0.2540 0.9837  0.0462 0.7232 0.2103 0.9798  0.0644 0.6908 0.2198 0.9750  0.0383 0.6688 0.2457 0.9528  0.0680 0.6753 0.2333 0.9766  End Members 0.0491^0.0485 0.6862^0.6770 0.2367^0.2570 0.9720^0.9825  0.0567 0.6889 0.2306 0.9761  0.0312 0.6112 0.3263 0.9687  0.0444 0.7002 0.2190 0.9636  0.0364 0.6593 0.2620 0.9577  0.0457 0.6706 0.2414 0.9577  Others  0.0059  0.0089  0.0077  0.0080  0.0065  0.0050^0.0056  0.0100  0.0044  0.0087  0.0072  0.0075  Fe 2 +  113  MM-39A MM-39A MM-39A MM-39A MM-39A MM-39 1^2^3^4^5^average  Table A4-2 Electron microprobe analyses of PLA amphiboles.  Sample Oxides  MM-293 1  MM-293 2  MM-299 1  MM-299 2  MM-29 average  MM-265 1  MM-265 2  MM-26 average  Si02 Al203 TiO2 FeO MgO MnO Na20 K2O CaO F Cl  47.80 5.99 0.97 12.16 16.08 0.81 1.41 0.30 11.41 0.42 0.13  49.86 5.21 0.90 12.16 16.24 0.86 1.35 0.22 11.43 0.32 0.06  41.39 12.96 3.05 8.70 15.94 0.06 2.91 0.65 12.33 0.32 0.02  40.49 13.23 2.98 8.58 15.88 0.74 2.89 0.65 12.32 0.45 0.02  44.88 9.35 1.98 10.40 16.03 0.62 2.14 0.46 11.87 0.38 0.06  49.95 5.33 1.02 11.45 17.00 0.82 1.30 0.25 11.28 0.34 0.08  48.94 5.93 0.97 11.46 16.54 0.78 1.43 0.32 11.28 0.46 0.13  49.44 5.63 1.00 11.46 16.77 0.80 1.36 0.28 11.28 0.40 0.10  Total  97.47  98.59  98.34  97.58  97.99  98.88  98.25  98.66  Si4 +  6.9192 1.0213 7.9405  Structural Formula (Based on 24 (0, OH, F, and Cl)) 7.2721 6.0533 5.8783 7.2746 7.0991 0.7279 1.9467 2.1217 0.7254 0.9009 8.0000 8.0000 8.0000 8.0000 8.0000  7.0991 0.9009 8.0000  3.4693 1.4717 0.0991 0.0000 0.1059 1.7700 0.0840 7.0000  3.5300 1.4828 0.1064 0.1679 0.0986 1.7866 0.0000 7.1723  3.4748 1.0635 0.0075 0.2865 0.3359 1.9322 0.0000 7.1003  3.4369 1.0424 0.0091 0.1426 0.3256 1.9163 0.1272 7.0000  3.6913 1.3945 0.1016 0.1899 0.1113 1.7605 0.0000 7.2491  3.5754 1.3909 0.0961 0.1135 0.1063 1.7532 0.0000 7.0353  3.5754 1.3909 0.0961 0.1135 0.1063 1.7532 0.0000 7.0353  Na+(A) K+  0.3117 0.0552 0.3669  0.3812 0.0399 0.4211  0.8264 0.1220 0.9484  0.6875 0.1207 0.8082  0.3666 0.0466 0.4131  0.4019 0.0591 0.4609  0.4019 0.0591 0.4609  OHFCl -  2.4421 0.1364 0.0265 2.6050  1.3650 0.1029 0.0121 1.4800  1.6215 0.1050 0.0032 1.7297  2.3455 0.1465 0.0037 2.4956  1.1547 0.1096 0.0151 1.2794  1.6948 0.1500 0.0258 1.8707  1.6948 0.1500 0.0258 1.8707  0  21.3950  22.5200  22.7297  21.5044  22.7206  22.1293  22.1293  A1 3-1- (IV) mg2+ Fe 2 + Mn2 + A1 3 +(VI) Ti4 + Ca2 + Na+(C)  114  Table A4-3 Electron microprobe analyses of PLA pyroxenes.  Sample MM-395 MM-395 MM- MM-39 MM-291 MM-291 MM-29 MM-407 MM-407 MM-407 MM-40 395 Oxides^2^3^4^average 1^2^average 1^2^3^average SiO2^54.27 54.50 54.09 54.28 52.07 50.16^51.12 54.20 54.48 52.80 53.82 Al 2 0 3^2.66^1.55^0.92^1.71^2.23^2.95^2.59^1.65^1.24^3.22^2.04 FeO^12.61^13.54^15.06^13.74 5.79^5.58^5.68^13.82^13.98^13.16^13.66 MgO^29.18^28.74 28.21 28.71 16.36^15.87^16.11 28.28^28.85^28.86^28.66 MnO^0.32^0.41^0.39^0.37^0.14^0.14^0.14 0.43^0.37^0.34^0.38 TiO 2^0.20^0.18^0.20^0.19 0.49^0.572 0.53 0.20^0.06^0.22^0.16 Na2 O^0.07^0.02^0.03^0.04 0.31^0.34^0.32 0.02^0.03^0.02^0.02 CaO^1.23^1.33^1.38^1.31^22.11^22.24^22.17 1.23^1.24^1.25^1.24 Total^100.54 100.26 100.27 100.36 99.50 97.84 ^98.67 99.83^100.26 99.88^99.99 Structural Formula (Based on 6 Oxygens) Si4 +^1.9240 1.9469 1.9482^1.9268 1.8931^1.9469 1.9502 1.8935 A1 3 ±(IV) 0.0760 0.0531 0.0392^0.0732 0.1069^0.0531 0.0498 0.1065 2.0000 2.0000 1.9875^2.0000 2.0000^2.0000 2.0000 2.0000 Ca2 +^0.0467 0.0508 0.0532^0.8765 0.8992^0.0474 0.0475 0.0482 Mg 2 +^1.5419 1.5303 1.5144^0.9020 0.8928^1.5141 1.5395 1.5428 Fe 2 +^0.3740 0.4047 0.4536^0.1793 0.1760^0.4153 0.4198 0.3946 M3+ (VI) 0.0351 0.0120 0.0000 ^0.0240 0.0242^0.0169 0.0027 0.0298 Mn 2 +^0.0095 0.0123 0.0119^0.0045 0.0043^0.0131 0.0112 0.0103 Na2 ÷^0.0050 0.0012 0.0023^0.0220 0.0252^0.0012 0.0020 0.0015 Ti4 +^0.0054 0.0049 0.0053^0.0137 0.0162^0.0054 0.0015 0.0060 2.0176 2.0163 2.0406^2.0220 2.0378^2.0133 2.0231 2.0332 End Members  Mg^0.7710 0.7652 0.7572^0.4510 0.4464^0.7570 0.7698 0.7714 Fe^0.1870 0.2023 0.2268^0.0896 0.0880^0.2077 0.2093 0.1973 Ca^0.0233 0.0254 0.0266^0.4383 0.4496^0.0237 0.0238 0.0241 0.9813 0.9929 1.0106^0.9789 0.9840^0.9884 1.0029 0.9928 Others^0.0275 0.0153 0.0097  ^  0.0321 0.0349^0.0183 0.0087 0.0238  115  Table A4-4 Original electron microprobe analyses of PLA biotites.  Sample Oxides  MM-265 core  MM-265 core  MM-264 rim  MM-264 rim  MM-264 core  MM-264 core  MM-26 average  Si02 Al203 FeO MgO MnO TiO2 K20 CaO Na20 F Cl  36.87 13.47 16.76 13.69 0.38 4.16 8.80 0.03 0.54 0.53 0.13  36.96 13.41 17.65 13.81 0.40 4.21 9.07 0.02 0.53 0.40 0.25  37.44 13.33 18.02 14.02 0.38 4.32 9.05 0.01 0.56 0.52 0.16  37.04 13.39 17.16 14.07 0.31 4.36 8.83 0.00 0.57 0.68 0.18  37.30 13.60 17.48 13.85 0.34 4.33 8.62 0.07 0.61 0.32 0.18  37.32 13.37 17.15 14.03 0.40 4.23 8.68 0.03 0.58 0.38 0.14  37.16 13.43 17.37 13.91 0.37 4.27 8.84 0.03 0.56 0.47 0.17  Total  95.37  96.62  97.80  96.59  96.69  96.30  96.56  Si 4 + Al 3 + (IV)  5.4779 2.3594 7.8373  A0+(V1) Mg2 + Fe 2 + Ti 4 + NIn 2-1-  0.0000 3.0334 2.0819 0.4650 0.0479 5.6283  0.0074 3.1234 2.2395 0.4808 0.0513 5.9024  0.1640 3.2104 2.3152 0.4988 0.0491 6.2375  0.0000 3.1619 2.1634 0.4941 0.0401 5.8595  0.0565 3.1195 2.2082 0.4922 0.0436 5.9200  0.0000 3.1438 2.1548 0.4782 0.0503 5.8266  K+ Na+ Ca2 +  1.6689 0.1545 0.0053 1.8288  1.7548 0.1570 0.0030 1.9148  1.7723 0.1674 0.0018 1.9416  1.6972 0.1670 0.0000 1.8643  1.6621 0.1785 0.0120 1.8526  1.6644 0.1704 0.0045 1.8393  OH FCl-  4.5878 0.1739 0.0267 4.7883  3.3286 0.1368 0.0525 3.5178  2.2524 0.1763 0.0336 2.4623  3.4327 0.2265 0.0381 3.6972  3.3325 0.1070 0.0366 3.4761  3.7084 0.1259 0.0290 3.8634  0  19.2117  20.4822  21.5377  20.3028  20.5239  20.1366  Structural Formula (Based on 24 (0, OH, F, and Cl)) 5.6087 5.7512 5.5837 5.6353 2.3913 2.2488 2.3794 2.3647 8.0000 8.0000 7.9631 8.0000  116  5.6075 2.3679 7.9754  Table A4.4 continued.  MM-294 MM-294 29 rim^rim^average  Sample Oxides  MM-297 core  MM-297 core  MM-297 rim  MM-297 rim  MM-294 MM-294 core core  SiO2 Al2°3 FeO MgO MnO TiO2 K2O CaO Na2O F Cl  37.31 13.23 16.45 14.50 0.32 4.30 8.58 0.07 0.54 0.73 0.20  36.76 13.34 16.07 14.48 0.30 4.39 8.25 0.11 0.53 0.72 0.20  37.15 13.43 16.33 14.55 0.35 4.21 8.83 0.10 0.60 0.60 0.20  36.90 13.42 15.74 14.71 0.29 4.34 8.65 0.04 0.70 0.53 0.23  36.95 13.29 17.18 13.81 0.28 4.34 8.55 0.09 0.52 0.48 0.19  37.13 13.57 17.58 13.75 0.33 4.232 8.64 0.04 0.51 0.64 0.18  36.80 13.32 17.11 13.72 0.38 4.27 8.59 0.08 0.50 0.62 0.15  36.76 13.32 17.27 14.01 0.39 4.31 8.38 0.01 0.49 0.64 0.18  36.97 13.36 16.72 14.19 0.33 4.3 8.56 0.07 0.55 0.62 0.19  Total  96.23  95.16  96.35  95.56  95.69  96.60  95.54  95.77  95.86  Si 4 + Al 3 +(IV)  5.5742 2.3302 7.9044  Structural Formula (Based on 24 (0, OH, F, and Cl)) 5.6000 5.4800 5.4162 5.5687 5.4685 5.5116 2.3157 2.3719 2.3450 2.3370 2.4000 2.3373 7.7319 7.9407 7.8135 7.8486 8.0000 7.8173  5.4865 2.3419 7.8284  Al 3 +(VI) Mg2 + Fe2 + 114+  0.0000 3.2288 2.0556 0.4827 0.0412 5.8053  0.0000 3.1806 1.9805 0.4868 0.0375 5.6854  0.0000 3.2507 2.0470 0.4744 0.0448 5.8186  0.0000 3.2505 1.9514 0.4839 0.0363 5.7221  0.0000 3.0170 2.1429 0.4873 0.0352 5.7363  0.0114 3.0914 2.2176 0.4799 0.0426 5.8428  0.0000 3.0459 2.1311 0.4778 0.0485 5.7034  0.0000 3.1154 2.1553 0.4839 0.0496 5.8043  K+ Na+ Ca2+  1.6343 0.1574 0.0114 1.8030  1.5508 0.1525 0.0177 1.7210  1.6892 0.1747 0.0155 1.8795  1.6354 0.2018 0.0069 1.8441  1.6275 0.1513 0.0139 1.7926  1.6625 0.1482 0.0070 1.8177  1.6311 0.1452 0.0123 1.7886  1.5960 0.1428 0.0018 1.7406  OH F  3.7581 0.2410 0.0408 3.0398  4.7513 0.2358 0.0400 5.0270  3.6455 0.2015 0.0409 3.8879  4.3821 0.1744 0.0472 4.6037  4.2949 0.1604 0.0394 4.4947  3.4162 0.2146 0.0368 3.6676  4.4296 0.2076 0.0310 4.6679  4.2077 0.2124 0.0379 4.4581  19.9602  18.9730  20.1121  19.3963  19.5053  20.3324  19.3321  19.5419  Mn 2 +  Ct  0  117  Table A4.4 continued.  Sample Oxides  MM-309 1  MM-309 MM-306 2^1  MM-306 MM-30 2^average  MM-392 core  MM-392 MM-392 rim^rim  MM-39 average  Si02 Al203 FeO MgO MnO TiO2 K20 CaO Na20 F Cl  36.29 12.49 15.99 13.73 0.37 4.05 7.93 0.17 0.59 1.09 0.16  36.02 13.02 16.44 13.74 0.37 4.18 7.90 0.29 0.59 0.50 0.15  36.74 13.60 16.94 14.01 0.36 4.22 8.79 0.04 0.71 0.36 0.24  37.24 13.73 17.03 13.53 0.48 4.40 8.67 0.05 0.70 0.48 0.19  36.57 13.21 16.60 13.75 0.40 4.21 8.32 0.14 0.65 0.61 0.19  36.34 14.19 16.85 13.96 0.25 4.32 8.68 0.03 0.58 1.38 0.16  36.62 13.20 17.14 13.87 0.51 4.27 9.02 0.00 0.64 0.98 0.16  37.13 13.36 17.16 14.06 0.29 4.35 8.94 0.01 0.58 0.98 0.14  36.69 13.58 17.05 13.96 0.35 4.31 8.88 0.01 0.60 1.11 0.15  Total  92.85  93.30  96.01  99.50  95.39  96.73  96.40  96.99  96.71  Si4 + A1 3 +(IV)  5.2084 2.1128 7.3212  Structural Formula (Based on 24 (0, OH, F, and Cl)) 5.2148 5.5111 5.6046 5.4610 5.5197 2.2215 2.4037 2.3954 2.5128 2.3456 7.4363 7.9149 8.0000 7.9738 7.8654  5.6166 2.3828 7.9994  Ti 4 + Mn 2 +  0.0000 2.9372 1.9187 0.4368 0.0452 5.3379  0.0000 2.9653 1.9904 0.4547 0.0454 5.4558  0.0000 3.1314 2.1246 0.4764 0.0454 5.7779  0.0395 3.0352 2.1437 0.4980 0.0612 5.7776  0.0000 3.1278 2.1182 0.4883 0.0321 5.7663  0.0000 3.1155 2.1612 0.4845 0.0648 5.8259  0.0000 3.1693 2.1714 0.4947 0.0374 5.8727  K+ Na+ Ca2+  1.4513 0.1636 0.0265 1.6413  1.4597 0.1657 0.0456 1.6710  1.6817 0.2078 0.0065 1.8960  1.6637 0.2030 0.0081 1.8748  1.6637 0.1677 0.0045 1.8359  1.7341 0.1842 0.0000 1.9213  1.7246 0.1692 0.0011 1.8950  FCl -  6.8534 0.3486 0.0319 7.2339  6.5585 0.1625 0.0297 6.7507  3.9879 0.1214 0.0503 4.1597  3.5168 0.1618 0.0398 3.7184  3.2789 0.4611 0.0329 3.7728  3.6125 0.3299 0.0335 3.9785  3.0390 0.3301 0.0287 3.3978  0  16.7661  17.2493  19.8403  20.2816  20.2272  20.0241  20.6022  A1 3 +(VI) Mg 2 + Fe2+  OH'  118  Table A4.4 continued.  Sample Oxides  MM-406 rim  MM-406 rim  MM-406 core  MM-406 core  MM-405 core  MM-405 core  MM-40 average  SiO2 Al203 FeO MgO MnO TiO2 K2O CaO Na2O F Cl  43.07 10.97 11.83 14.07 0.38 2.52 0.26 11.46 2.37 0.17 0.01  43.88 10.49 11.31 14.64 0.22 2.36 0.28 11.92 2.39 0.23 0.05  41.55 13.09 12.72 13.22 0.24 2.36 0.51 11.76 2.24 0.32 0.06  43.28 11.94 11.15 14.15 0.33 2.30 0.57 11.26 2.02 0.13 0.09  36.23 13.55 16.62 13.92 0.38 4.43 8.76 0.06 0.64 0.75 0.15  36.98 13.29 16.77 14.19 0.35 4.27 8.80 0.03 0.68 0.46 0.16  40.83 12.22 13.40 14.03 0.32 3.04 3.20 7.75 1.73 0.34 0.09  Total  97.10  97.77  98.08  97.22  95.48  95.99  96.94  Si4 + A1 3-1- (IV)  6.2533 1.7467 8.0000  A1 31- (VI) Mgt -1Fe2 + Ti 4 + Mn2 +  0.1300 3.0440 1.4362 0.2757 0.0472 4.9331  0.2004 3.1826 1.3787 0.2583 0.0273 5.0472  0.4082 2.9073 1.5699 0.2614 0.0294 5.1761  0.2983 3.0538 1.3493 0.2505 0.0402 4.9891  0.0000 3.0876 2.0686 0.4957 0.0474 5.6992  0.0000 3.1666 2.0991 0.4810 0.0447 5.7914  K+ NalCa2+  0.0477 0.6675 1.7822 2.4975  0.0513 0.6768 1.8623 2.5904  0.0957 0.6426 1.8646 2.6029  0.1056 0.5680 1.7458 2.4194  1.6634 0.1863 0.0091 1.8588  1.6815 0.1984 0.0047 1.8846  OH FCl -  2.8042 0.0546 0.0020 2.8608  2.1642 0.0756 0.0109 2.2507  1.8914 0.1046 0.0120 2.0080  2.6736 0.0412 0.0186 2.7334  4.4892 0.2476 0.0310 4.7679  4.0076 0.1544 0.0348 4.1968  0  21.1392  21.7493  21.9920  21.2666  19.2321  19.8032  -  Structural Formula (Based on 24 (0, OH, F, and Cl)) 6.3973 6.1316 6.2623 5.3919 1.6027 1.8684 1.7377 2.3771 8.0000 8.0000 8.0000 7.7690  119  5.5376 2.3448 7.8823  Input Data File In.Dat An Assimilation Example from Literature. 0.005 9 3 Oxide Parent Daughter Plag 57.76 58.65 53.35 Si02 0.89 0.77 0.04 TiO2 Al203 17.18 17.42 29.51 6.46 6.06 0.93 FeO 5.61 4.00 0.10 MgO 6.90 6.41 CaO 11.80 3.69 3.96 4.42 Na2O K20 1.22 1.49 0.23 0.29 0.30 0.00 P205 0.90  Opx 53.86 0.00 2.46 13.64 28.53 1.48 0.00 0.00 0.00  Xenolith 72.01 0.20 14.44 1.78 0.89 2.61 3.96 3.00 0.13  Relative_Error 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0  Description of Data File IN.DAT Line 1:^(Characters) The title of the project. Line 2:^(Float) The interval of amounts of minerals to map the x 2 and R2 distributions around the optimal solution. Line 3:^(Integers) The numbers of the oxide components (n) and mineral phases (m). Line 4:^(Character strings) Labels (such as "oxide"), the name of the parent rock, the name of derivative rock, and the names of (m) mineral phases. Note: Each character string must not contain any spaces. Line 5 - (4+n): (One character string, and m+3 floats) The name of oxide component, the compositions of the parent and derivative rocks, the compositions of the mineral phases, and the relative analytical error for each oxide. Line (5+n):^(Float) The confidence level used to calculate the confidence region.  120  Example of Ouput of Program MASSBAL.C. Filename = Out.dat.  ************ ******************* **** An Assimilation Example from Literature ***********************************  INPUT PARAMETRS number of oxide = 9,^number of mineral_phase = 3,^degrees_of freedom = 6 Oxide SiO2 TiO2 Al203 FeO MgO CaO Na2O K2O P2O5  Parent 57.76 0.89 17.18 6.46 5.61 6.90 3.69 1.22 0.29  Daughter 58.65 0.77 17.42 6.06 4.00 6.41 3.96 1.49 0.30  Plag 53.35 0.04 29.51 0.93 0.10 1.80 4.42 0.23 0.00  Opx 53.86 0.00 2.46 13.64 28.53 1.48 0.00 0.00 0.00  Xenolith 72.01 0.20 14.44 1.78 0.89 2.61 3.96 3.00 0.13  Relative Error 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00  The Parameters of the Compositional Difference Matrix # 1 2 3 4 5 6 7 8 9  oxide SiO2 TiO2 Al203 FeO MgO CaO Na2O K2O P2O5  ml-d (-5.30 (-0.73 (12.09 (-5.13 (-3.90 (5.39 (0.46 (-1.26 (-0.30  m2-d -4.79 -0.77 -14.96 7.58 24.53 -4.93 -3.96 -1.49 -0.30  m3-d 13.36) -0.57) -2.98) -4.28) -3.11) -3.80) 0.00) 1.51) -0.17)  =  = = = =  = = = =  p-d -0.89 0.12 -0.24 0.40 1.61 0.49 -0.27 -0.27 -0.01  The Parameters of the Weighted [to 1/standard_deviation] Matrix # 1 2 3 4 5 6 7 8 9  Oxide SiO2 TiO2 Al203 FeO MgO CaO Na2O K2O P2O5  ml-d (-9.176 (-82.022 (70.373 (-79.412 (-69.519 (78.116 (12.466 (-103.279 (-103.448  m2-d -8.293 -86.517 -87.078 117.337 437.255 -71.449 -107.317 -122.131 -103.448  m3-d 23.130) -64.045) -17.346) -66.254) -55.437) -55.072) 0.000) 123.770) -58.621)  121  = =  =  = =  = =  = =  pd -1.541 13.483 -1.397 6.192 28.699 7.101 -7.317 -22.131 -3.448 -  Standard deviation 0.578 0.009 0.172 0.065 0.056 0.069 0.037 0.012 0.003  The x2 Solution.  The Parameters of the Matrix U from SVDCMP 0.020 0.133 0.175 -0.255 -0.860 0.138 0.207 0.243 0.162  (^  (^  (^  (  (  (^  (^  (^  (^  0.052 0.388 -0.286 0.289 0.110 -0.339 -0.010 0.557 0.494  -0.122^ 0.457^ 0.144^ 0.311^ -0.000^ 0.345^ 0.079^ -0.579^ 0.446^  )  )  )  )  )  )  )  )  )  The Parameters of the Diagonal Matrix W 518.23968506 219.75834656 173.74630737 The Parameters of the Matrix V 0.102 -0.986 0.128  (^  (  (^  -0.993 -0.094 0.067  -0.054 -0.134 -0.989  )  )  )  122  The Minimized x2 Solution (out of the total 1.0 gram) mineral_phase weight ^ 0.0137051 grams plag ^ 0.0455283 grams Opx Xeno^-0.1302265 grams Weighted Residual on Each Oxide Equation = RHS - LHS  # 1 2 3 4 5 6 7 8 9  oxide Si02 TiO2 Al203 FeO MgO CaO Na20 K20 P205  residual 1.974620 10.205893 -0.655798 -6.689886 2.524708 2.111935 -2.601957 0.962923 -4.954653  Chi_Squared Value = 196.324387 Residual on Each Oxide Equation = RHS - LHS (using unweighted compositions)  # 1 2 3 4 5 6 7 8 9  Oxide Si02 TiO2 Al203 FeO MgO CaO Na2O K20 P205  residual 1.140541 0.090832 -0.112666 -0.432167 0.141636 0.145724 -0.096012 0.011748 -0.014368  The R2 Associated with the x 2 Solution= 1.559404  Del-Chi-Squared Mapped against Amounts of Phases  Plag vs Opx (Xenolith = -0.1302 gram Plag \ Opx -0.006 -0.001 0.004 0.009 0.014 0.019 0.024 0.029 0.034  0.026 107.33 102.94 101.08 101.74 104.93 110.65 118.90 129.67 142.97  0.031 65.87 60.37 57.40 56.95 59.03 63.63 70.76 80.42 92.60  0.036 37.54 30.92 26.83 25.27 26.23 29.72 35.74 44.29 55.36  0.041 22.32 14.59 9.38 6.71 6.56 8.94 13.84 21.27 31.22  0.046 20.21 11.37 5.05 1.26 0.00 1.26 5.05 11.37 20.21  0.051 31.22 21.27 13.84 8.94 6.56 6.71 9.38 14.59 22.31  0.056 55.36 44.29 35.74 29.72 26.23 25.27 26.83 30.92 37.54  0.061 92.60 80.42 70.76 63.63 59.03 56.95 57.40 60.37 65.87  0.066 142.97 129.67 118.90 110.65 104.93 101.74 101.08 102.94 107.33  Plag vs Xenolith (Opx = 0.0455 gram ) Plag Xeno -0.006 -0.001 0.004 0.009 0.014 0.019 0.024 0.029 0.034  -0.150 35.42 26.19 19.49 15.32 13.68 14.56 17.96 23.90 32.36  -0.145 29.05 19.92 13.32 9.24 7.69 8.67 12.17 18.20 26.76  -0.140 24.40 15.36 8.85 4.87 3.42 4.49 8.09 14.21 22.86  -0.135 21.45 12.51 6.10 2.21 0.85 2.02 5.72 11.94 20.68  -0.130 20.21 11.37 5.05 1.26 0.00 1.26 5.05 11.37 20.21  -0.125 20.68 11.94 5.72 2.02 0.85 2.21 6.10 12.51 21.45  -0.120 22.86 14.21 8.09 4.49 3.42 4.87 8.85 15.36 24.40  -0.115 26.76 18.20 12.17 8.67 7.69 9.24 13.32 19.92 29.05  -0.110 32.36 23.90 17.96 14.56 13.68 15.32 19.49 26.19 35.42  -0.115 130.79 80.34 43.01 18.79 7.69 9.71 24.84 53.10 94.46  -0.110 142.83 90.86 52.02 26.29 13.68 14.18 27.80 54.54 94.39  Opx vs Xenolith (Plag = 0.0137 gram ) Opx \ Xeno 0.026 0.031 0.036 0.041 0.046 0.051 0.056 0.061 0.066  -0.150 94.39 54.54 27.80 14.18 13.68 26.29 52.02 90.86 142.83  -0.145 94.46 53.10 24.84 9.71 7.69 18.79 43.01 80.34 130.79  -0.140 96.24 53.36 23.60 6.95 3.42 13.00 35.71 71.53 120.46  -0.135 99.73 55.34 24.06 5.90 0.85 8.93 30.12 64.42 111.84  -0.130 104.93 59.03 26.23 6.56 0.00 6.56 26.23 59.03 104.93  -0.125 111.84 64.42 30.12 8.93 0.85 5.90 24.06 55.34 99.73  -0.120 120.46 71.53 35.71 13.00 3.42 6.95 23.60 53.36 96.24  ****************************************************************  124  The R 2 Solution.  The Paramenters of the Matrix U from SVDCMP  (^0.115 (^0.013 (^0.546 ( -0.281 ( -0.748 (^0.188 (^0.118 (^0.033 (^0.005  0.304 -0.181 0.405 -0.606 0.551 0.051 -0.155 -0.096 -0.068  -0.822^) -0.002^) 0.405^) 0.095^) 0.204^) 0.308^) -0.015^) -0.120^) -0.004^)  The Parameters of the Diagonal Matrix W 32.80237579 6.23699284 17.70904922 The Parameters of the Matrix V  (^0.346 ( -0.934 (^0.084  0.758 0.332 0.561  0.552 0.131 -0.823  ) ) )  The Minimized R 2 Solution (out of the total 1.0 gram) Mineral_phase weight ^ Plag 0.0636599 grams ^ Opx 0.0699627 grams Xeno^-0.0276482 grams Residual on Each Oxide Equation = RHS - LHS  # 1 2 4 5 6 7 8 9  Oxide Si02 TiO2 FeO MgO CaO Na20 K20 P205  residual 0.151895 0.204584 0.077924 0.056103 0.386727 -0.022231 -0.043795 0.025387  The Least Residual_Squared = 0.228821  125  Del-Residual-Squared Mapped against Amounts of Phases  Plag vs Opx (Xenolith = -0.0276 grams ) Plag\Opx 0.044 0.049 0.054 0.059 0.064 0.069 0.074 0.079 0.084  0.050 0.226 0.246 0.278 0.323 0.380 0.449 0.531 0.625 0.731  0.055 0.123 0.127 0.143 0.172 0.213 0.267 0.333 0.411 0.502  0.060 0.067 0.056 0.056 0.069 0.095 0.133 0.183 0.245 0.320  0.065 0.059 0.032 0.017 0.014 0.024 0.046 0.080 0.127 0.186  0.070 0.099 0.056 0.025 0.006 0.000 0.006 0.025 0.056 0.099  0.075 0.186 0.127 0.080 0.046 0.024 0.014 0.017 0.032 0.059  0.080 0.320 0.245 0.183 0.133 0.095 0.069 0.056 0.056 0.067  0.085 0.502 0.411 0.333 0.267 0.213 0.172 0.143 0.127 0.123  0.090 0.731 0.625 0.53 0.449 0.380 0.323 0.278 0.246 0.226  Plag vs Xenolith (Opx = 0.0700 grams ) Plag\Xeno 0.044 0.049 0.054 0.059 0.064 0.069 0.074 0.079 0.084  -0.048 0.116 0.092 0.080 0.080 0.093 0.118 0.156 0.205 0.268  -0.043 0.094 0.065 0.049 0.044 0.052 0.073 0.105 0.151 0.208  -0.038 0.084 0.050 0.029 0.020 0.023 0.039 0.067 0.107 0.160  -0.033 0.086 0.047 0.021 0.007 0.006 0.017 0.040 0.076 0.124  -0.028 0.099 0.056 0.025 0.006 0.000 0.006 0.025 0.056 0.099  -0.023 0.124 0.076 0.040 0.017 0.006 0.007 0.021 0.047 0.086  -0.018 0.160 0.107 0.067 0.039 0.023 0.020 0.029 0.050 0.084  -0.013 0.208 0.151 0.105 0.073 0.052 0.044 0.049 0.065 0.094  -0.008 0.268 0.205 0.156 0.118 0.093 0.080 0.080 0.092 0.116  Opx vs Xenolith (Plag = 0.0637 grams ) Opx\Xeno 0.050 0.055 0.060 0.065 0.070 0.075 0.080 0.085 0.090  -0.048 0.384 0.240 0.143 0.094 0.093 0.139 0.232 0.373 0.561  -0.043 0.365 0.216 0.114 0.059 0.052 0.093 0.181 0.316 0.499  -0.038 0.358 0.203 0.096 0.036 0.023 0.058 0.140 0.270 0.447  -0.033 0.363 0.203 0.090 0.024 0.006 0.035 0.112 0.236 0.408  -0.028 0.380 0.213 0.095 0.024 0.000 0.024 0.095 0.213 0.380  126  -0.023 0.408 0.236 0.112 0.035 0.006 0.024 0.090 0.203 0.363  -0.018 0.447 0.270 0.140 0.058 0.023 0.036 0.096 0.203 0.358  -0.013 0.499 0.316 0.181 0.093 0.052 0.059 0.114 0.216 0.365  -0.008 0.561 0.373 0.232 0.139 0.093 0.094 0.143 0.240 0.384  ************************************************************************** Del Chi_Sq Value = 10.644501 less than which the Probability of the Chi_Sq Variable Is 0.90000 ************************************************************************** The Orientations and Lengths of Semi-Axes of the x2 Solution Confidence Region (0.102, -0.986, 0.128) (-0.993, -0.094, 0.067) (-0.054, -0.134, -0.989)  semi-axis length = 0.0063 semi-axis length = 0.0148 semi-axis length = 0.0188  127  /*^ /*^ /*^  PROGRAM MASSBAL.0 BY: PEIWEN KE May, 1992  /* This program is written to calculate the amount of minerals fractionating from magma A to magma B. The principles involved are material balance, and the Chi-squared minimization fitting technique. */ #include <stdio.h> #include <math.h> #include "ganunln.c" #include "gammq.c" #include "gcf.c" #include "gser.c" #include "malloc.h" #include "nrutil.c" #include "nrutil.h" #define TOL 1.0e-5 /* THE MAIN PROGRAM MASSBAL.0  ^  */  main(argc, argv) int argc; char **argv; {  int i,j, k, il,j1, kl, I; int number_of oxide; int number_of mineral_phase; int degrees_of freedom; int number_of characters_in_title; int string_length(); char title[81]; char name_of component[10]; char name_of oxide[12][12]; char name_of_parent[20]; char name_of daughter[20]; char name_of mineral_phase[20][8]; char name_of error[20]; float step, alpha, r_2, df; float wmax, thresh, sum, sum_chi_to_r, chisq, residual_square, residual_chi_to_r; float mole_pt_of_p[12]; float mole_pt_of d[12]; float mole_pt_of m[12][8]; float relative_error[12]; float std_deviation[12]; float diff of_p_and_d[12]; float diff of m_and_d[12][8]; float diff of m_and_d_save[12][8]; float diff of_p_and_d_std[12]; float diff of m_and_d_std[12][8]; float diff of m_and_d_std_save[12][8]; 128  float weight_of mineral_phase[8]; float weight_of mineral_phase_r[8]; float residual_std[12]; float residual_std_chi_to_r[12]; float residual_412]; float w[8]; float v[8][8]; float weight_of mineral[8][16]; float chisq_array[16][16][8]; float chisq_diff array[16][16][8]; float one_phase_chisq_diff array[16]; float one_phase_residual_square_diff array[16]; float w_r[8]; float v_r[8][8]; float weight_of mineral_r[8][16]; float residual_square_array[16][16][8]; float residual_square_diff array[16][16][8]; float p; float one_minus_p; float half degrees_of freedom; float left_ side_ of_ root; float nght side_of root; float fittedidelta_chisq; float left_side_of fitted_delta_chisq; float right_side_of fitted_delta_chisq; float motification_of theleft_and_right_sides; float delta_weight_of mineral_phase[8]; float delta; float unit_length_of axis[8]; float delta_chisq_length_of axis[8]; .  FILE *fp; FILE *fo; FILE *fo2; fp = fopen(argv[1], "r"); fo = fopen(argv[2], "w"); fo2 = fopen("plot.out", "w");  /* inputting data file */ /* outputting data file */ /* outputting data file for plotting only for three mineral phases*/  if (argc != 3) {  )  {  printf("please check the opening data file\n"); exit(0); else if ((fp = fopen(argv[1], "r")) = = NULL) printf("Can't open file %s\n", argv[1]); exit(0);  }  /* initiate the parameters used in calling the functions */ for (i=0; i< =numberof mineral_phase; i+ +) 129  {  w[i]=0.00; for (j =0; j < =number of mineral_phase; j + +) v[i][j]=0.00; —  }  for (i = 0; i< = number_of mineral_phase; i = i +1) weight_of mineral_phase[i] = 0.0; /* input the raw data */ fscanf (fp, "%[^ \n]", title); fscanf (fp, " %An", &step); fscanf (fp, "%d%d\n“, &number_of oxide, &number_of mineral_phase); fscanf (fp, " %s %s %s ", name_of component, name_of_parent, name_of_daughter); for (i=1; i< =number_of mineral_phase; i+ +) fscanf (fp, "%s", name_of mineral_phase[i]); fscanf (fp, "%s\n", name of error); number_of characters_in_title = string_length (title); degrees_of freedom = number_of_oxide - number_of_mineral_phase; half degrees_of freedom = degrees_of freedom/2.0; for (i =1; i< = 2 +number_of characters_intitle; i+ +) fprintf (fo, "*"); fprintf (fo, "\n"); fprintf (fo, " %s \n", title); for (i =1; i< =2 +number_of characters_in_title; i+ +) fprintf (fo, "*"); fprintf (fo, "\n\n\n"); fprintf (fo, "INPUT PARAMETRS\n"); fprintf (fo, " ^\n"); fprintf (fo, "number_of oxide = %d, number_of_mineral_phase = %d, degrees_of freedom = %d\n", number_of_oxide, number_of_mineral_phase, degrees_of freedom); fprintf (fo, " % 10s % 10s % 10s", name_of component, name_of_parent, name_of_daughter); for (i=1; i< =number_of mineral_phase; i+ +) fprintf (fo, " % 10s", name_of mineral_phase[i]); fprintf (fo, " %18s\n",name_of error); for (i=1; i< =number_of oxide; i+ +) {  fscanf (fp, "%s%f%f", name_of oxide[i], &mole_pt_of_p[i], &mole_pt_of d[i]); fprintf (fo, " % 10s %10.2f%10.2f", name_of oxide[i], mole_pt_of_p[i], mole_pt_of d[i]); diff of_p_and_d[i] = mole_pt_of_p[i] - mole_pt_of d[i]; for (j =1; j < =number_of mineral_phase; j+ +) {  fscanf (fp, "%f", &mole_pt_of m[i][j]); fprintf (fo, "%10.2f ", mole_pt_of m[i][j]); diff of m_and_d[i][j] = mole_pt_of m[i][j] - mole_pt_of d[i]; diff of mand_dsave[i][j] = diff of mandd[i][j]; }  fscanf (fp, "%f", &relative_error[i]); 130  ^  fprintf (fo, "%15.2f", relativeerror[i]); std_deviation[i] = 1.0; if (mole_pt_of_p[i] != 0.00) if (relative_error[i] != 0.00) std_deviation[i] = mole_pt_of_p[i]*relative_error[i]/100.0; diff of_p_and_d_std[i] = diff of_p_and_d[i]/std_deviation[i]; for (j = 1; j < =number_of mineral_phase; j + +) { diff of m_and_d_std[i][j] = diff of m_and_d[i][Wstd_deviation[i]; diff of m_and_d_std_save[i][j]=diff of m_and_d_std[i][j]; } fprintf (fo, " \n"); }  fscanf (fp, "%f", &p); one_minus_p = 1.0 - p; /* output the parameters of the compositional difference matrix */ fprintf (fo, "\nThe Parameters of the Compositional Difference Matrix\n"); fprintf (fo, " ^ \n"); fprintf (fo, "# oxide^"); for (i= 1; i< =number_of mineral_phase; i=i +1) fprintf(fo, "m%i-d ", i); fprintf (fo, " p-d \n"); for (i= 1; i< =number_of oxide; i=i +1) { fprintf (fo, "%i %6s ( ", i, name_of oxide[i]); for (j =1; j < =number_of mineral_phase; j + +) fprintf (fo, " %5.2f ", diff of m_and_d[i][j]); fprintf (fo, ") = %5.2f\n", diff of_p_and_d[i]); } /* output the parameters of the weighted [to 1/std_deviation] matrix A */ fprintf (fo, "\nThe Parameters of the weighted to [1/standard_deviation] Matrix\n"); fprintf (fo, " ^ \n"); fprintf (fo, "# oxide^"); for (i=1; i< =number_of mineral_phase; i+ +) fprintf(fo, "m%i-d^", i); fprintf (fo, "^p-d^standard deviation\n"); for (i= 1; i < =number_of oxide; i + +) {  fprintf (fo, "%i %6s ( ", i, name_of oxide[i]); for (j = 1; j < =number_of mineral_phase; j + +) fprintf (fo, "%7.3f ", diff of m_and_d_std[i][j]); fprintf (fo, ") = %7.3f^%7.3f\n", diff of_p_and_d_std[i], std_deviation[i]);  }  svdcmp(diff of m_and_d_std, number_of oxide, number_of mineral_phase, w, v); wmax =0.0; 131  for (j =1; j < = number_of mineral_phase; j+ +) if (w[j] > wmax ) wmax=w[j]; thresh = TOL *wtnax; for (j =1; j < =number of mineral_phase; j + +) if (w[j] < thresh) w[j]=0.0; fprintf (fo, "\n\nThe Chi-Squared Solution. \n"); fprintf (fo, " ^  \n");  1* output the parameters of the U matrix from svdcmp *1 fprintf (fo, "\nThe Parameters of the Matrix U from SVDCMP\n"); fprintf (fo, " ^ \n"); for (i=1; i< =number_of oxide; i= i + 1) {  fprintf (fo, "( "); for (j =1; j < =number of mineral_phase; j + +) fprintf (fo, " %7.3f ", diff of m_and_d_std[i][j]); fprintf (fo, ")\n"); }  1* output the parameters of the diagonal matrix W *1 fprintf (fo, "\nThe Parameters of the Diagonal Matrix W\n"); fprintf (fo, " ^ \n"); for (i=1; i< =numberof mineral_phase; i=i +1) fprintf (fo, " % 12.8f ", w[i]); fprintf (fo, "\n"); 1* output the parameters of the matrix V *1 fprintf (fo, "\nThe Parameters of the Matrix V \n"); fprintf (fo, " ^\n"); for (i =1; i< =numberof mineral_phase; i = i +1) {  fprintf (fo, "( "); for (j =1; j < =numberof mineral_phase; j=j + 1) fprintf (fo, " %7.3f ", v[i][j]); fprintf (fo, " )\n"); }  svbksb(diff of m_and_dstd, w, v, number_of oxide, number_of mineral_phase, diff of_p_and_d_std, weight_of mineral_phase); chisq = 0.0; residual_chi_tor = 0.0; for (i=1; i< =number_of oxide; i + +) {  sum = 0.0; sumchito_r = 0.0; for (j =1; j < =number_of mineral_phase; j + +){ sum + = diff of m_and_d_std_save[i][j]*weight_of mineral_phase[j]; sum_chi_to_r += diff of m_and_d_save[i][j]*weight_of mineral_phase[j];} residual_std[i] = diff of_p_and_d_std[i] - sum; 132  residual_std_chi_to_r[i] = diff of_p_and_d[i] - sum_chi_to_r; chisq + = residual_std[i]*residual_std[i]; residual_chi_to_r + = residual_std_chi_to_r[i]*residual_std_chi_to_r[i]; }  /* output of the solution of least chi-square */ fprintf (fo, " \nThe Minimized Chi-Squared Solution (out of the total 1.0 gram)\n"); fprintf (fo, " ^ \n"); fprintf (fo, "^mineral_phase^weight \n"); for (i= 1; i< =number_of mineral_phase; i=i + 1) fprintf(fo, " %20s % 10.7f grams \n", name_of mineral_phase[i], weight_of mineral_phase[i]); fprintf (fo, " \n"); /* output of residual on each oxide equation */ fprintf (fo, "Weighted Residual on Each Oxide Equation = RHS - LHS\n"); fprintf (fo, " ^ \n"); fprintf (fo, "# oxide^residual \n"); for (i= 1; i< =number_of oxide; i+ +) fprintf (fo, "%d %6s^%f \n", i, name_of oxide[i], residual_std[i]); fprintf (fo, " \n"); /* output the least chi_square */ fprintf (fo, "Chi_Squared Value = %fln", chisq ); fprintf (fo, "Residual on Each Oxide Equation = RHS - LHS (using unweighted compositions)\n"); fprintf (fo, " ^ \n"); fprintf (fo, "# oxide^residual \n"); for (i= 1; i< =number_of oxide; i+ +) fprintf (fo, "%d %6s^%f \n", i, name_of oxide[i], residual_std_chi_to_r[i]); fprintf (fo, "\n"); fprintf (fo, "The Residual-Squared Associated with the Chi_Squared solution= %fin", residual_chi_to_r ); for (i1=1; il < =number_of mineral_phase; il + +) for (i=1; i< =9; i++) weight_of mineral[il][i] = weight_of mineral_phase[il] + step*(i-5); /* output the distribution of the chi_square difference between the best solution and the chi_square around it */ fprintf (fo, " \nDel-Chi-Sqaured Mapped against Amounts of Phases\n"); fprintf (fo, " ^ \n"); 1=0; for (i= 1; i< =number_of mineral_phase; i + +) if (number_of mineral_phase= =1) {  fprintf (fo, " % lOs ", name_of mineral_phase[i]); for (j=1; j<=9; j++) 133  fprintf (fo, " %6.3f ", weight_of mineral[1][j]); fprintf (fo, "\n"); fprintf (fo, "delta chisq "); for (j=1; j<=9; j++){ one_phase_chisq_diff array[j] = 0.0; {  for (k =1; k< =number_of oxide; k+ +) onephase_chisq_diff array[j] + = ((diff ofp_and_d_std[k] diff of m_and_d_std_save[k][1]*weight_of mineral[1][j])*(diff of_p_and_d_std[k] diff of m_and_d_std_save[k][1]*weight_of mineral[1][j]) - chisq); —  }  fprintf (fo, " %6.3f ", one_phase_chisq_diff array[j]); } fprintf (fo, "\n"); } else {  for (j =i + 1; j < =numberof mineral_phase; j + +){ I+ =1; fprintf (fo, " \n %s vs %s (", name_of mineral_phase[i], name_of mineral_phase[j]); for (k =1; k< =numberof mineral_phase; k+ +){ if (k!=i) if (k! =j) fprintf (fo, "%s = %7.4f gram ", name of mineral_phase[k], weight_of mineralphase[k]); }  fprintf (fo, ")\n"); fprintf (fo, " %s\\ %s ", name_of mineral_phase[i], name_of mineral_phase[j]); for (i1=1; il< =9; il++) fprintf (fo, "%6.3f ", weight_of mineral[j][il]); fprintf (fo, "\n"); for (il =1; il < =9; il+ +) fprintf (fo, " %6.3f ", weight_of mineral[i][il]); for (j1=1; jl<=9; j1++){ chisq_array[il][jl][I]=0.0; for (k =1; k< =numberof oxide; k+ +) { sum = 0.0; for (kl =1; kl < =numberof mineral_phase; kl + +){ if (kl= =i) sum += diff of m_and_d_std_save[k][kl]*weightof mineral[kl][il]; else { if (kl= =j) sum += diff of m_and_d_std_save[k][k1]*weight_of mineral[kl][j1]; else sum + = diff of mand_d_stdsave[k][k1]*weight_of mineralphase[kl]; }  }  residual_std[k] = diff of_p_and_d_std[k] - sum; chisq_array[i1][j1][I] + = residual_std[k]*residual_std[k]; }  chisq_diff array[i1] [j 1] [I] = chisq_array[il][jl][I]-chisq; fprintf (fo, " %6. 2f ", chisq_diff array[i 1 ][j MID; }  134  fprintf (fo, " \n"); }  }  }  /* */  /* output the calculation results based on the residual_square */ svdcmp(diff of m_and_d, number_of oxide, number_of mineral_phase, w_r, v_r); wmax =0.0; for (j =1; j < =number_of mineral_phase; j+ +) if (w_r[j] > wmax ) wmax =w_r[j]; thresh = TOL *wmax; for (j =1; j < =number_of mineral_phase; j + +) if (w_r[j] < thresh) w_r[j]=0.0; fprintf (fo, nn ****************************************************************\ n\n ”) ; fprintf (fo, "The Residual_Squared Solution\n"); fprintf (fo, " ^ \n"); /* output the parameters of the U matrix from svdcmp from r-square */ fprintf (fo, " \nThe Paramenters of the Matrix U from SVDCMP\n"); fpiintf (fo, " ^ \n"); for (i = 1; i< =number_of oxide; i=i + 1) {  fprintf (fo, "( "); for (j =1; j < =number_of tnineral_phase; j + +) fprintf (fo, " %7.3f ", diff of m_and_d[i][j]); fprintf (fo, ")\n"); }  /* output the parameters of the diagonal matrix W based on r_square */ fprintf (fo, " \nThe Parameters of the Diagonal Matrix W\n"); ^ fprintf (fo, " \n"); for (i=1; i< =number_of mineral_phase; i=i +1) fprintf (fo, " %12.8f ", w_r[i]); fprintf (fo, " \n"); /* output the parameters of the matrix V based on r_square */ fprintf (fo, " \nThe Parameters of the Matrix V \n"); ^ fprintf (fo, " \n"); for (i = 1; i < =number_of mineral_phase; i=i + 1) {  fprintf (fo, "( "); for (j = 1; j < =number_of mineral_phase; j=j +1) fprintf (fo, " %7.3f ", v_r[i][j]); fprintf (fo, " )\n"); 135  }  svbksb(diff of m_and_d, w_r, v_r, number_of oxide, number_of mineral_phase, diff of_p_and_d, weight_of mineral_phase_r); residual_square = 0.0; for (i= 1; i< =number_of oxide; i+ +) {  sum = 0.0; for (j = 1; j < =number_of mineral_phase; j + +) sum + = diff of m_and_d_save[i][j]*weight_of mineral_phase_r[j]; residual_r[i] = diff of_p_and_d[i] - sum; residual_square + = residual_r[i]*residual_r[i]; }  /* output of the solution of least residual-square */ fprintf (fo, "\nThe Minimized Residual-Squared Solution (out of the total 1.0 gram)\n"); fprintf (fo, " ^ \n"); fprintf (fo, "^Mineral_phase^weight \n"); for (i = 1; i< =number_of mineral_phase; i= i + 1) fprintf(fo, " %20s^%10.7f grams \n", name_of mineral_phase[i], weight_of mineral_phase_r[i]); fprintf (fo, " \n"); /* output of residual on each oxide equation */ fprintf (fo, "Residual on Each Oxide Equation = RHS - LHS\n"); fprintf (fo, " ^ \n"); fprintf (fo, "# oxide^residual \n"); for (i= 1; i< =number_of oxide; i+ +) fprintf (fo, "%d %6s^%f \n", i, name_of oxide[i], residual_r[i]); fprintf (fo, " \n"); /* output the least residual_square */ fprintf (fo, "The Least Residual_Squares = %An", residual_square ); for (i1=1; il < =number_of mineral_phase; il + +) for (i=1; i< =9; i++) weight_of mineral_r[il][i] = weight_of mineral_phase_r[il] + step*(i-5); /* output the distribution of the residual_square difference between the best solution and the residual_square around it */ fprintf (fo, "\nDel-Residual-Squared Mapped against Amounts of Phases\n"); fprintf (fo, " ^ \n"); I=0; for (i = 1; i< =number_of mineral_phase; i+ +) if (number_of mineral_phase= =1) {  fprintf (fo, " % 10s ", name_of mineral_phase[i]); 136  for(j=1; j<=9; j++) fprintf (fo, " %6.3f ", weightof mineral_r[l][j]); fprintf (fo, "\n"); fprintf (fo, "delta-resisq "); for(j=1; j<=9; j++){ one_phase_residual_square_diff array[j] = 0.0; {  for (k= 1; k< =number_of oxide; k+ +) one_phase_residual_square_diff array[j] + = ((diff of_p_and_d[k] diff of m_and_d_save[k][1]*weight_of mineral_r[l][j])*(diff of_p_and_d[k] diff of m_and_d_save[k][1]*weight_of mineral_r[l][j]) - residual_square); fprintf (fo, " %6.3f ", one_phase_residual_square_diff array[j]); }  fprintf (fo, "\n"); }  else {  for (j=i +1; j< =number_of mineral_phase; j + +){ I+ =1; fprintf (fo, " \n %s vs %s (" , name_of mineral_phase[i], name_of mineral_phase[j]); for (k=1; k< =numberof mineral_phase; k+ +){ if (k!=i) if (k!=j) fprintf (fo, "%s = %7.4f grams ", name_of mineral_phase[k], weight_of mineral_phase_r[k]);} fprintf (fo, ")\n"); fprintf (fo, " %s\\ %s ",name_of mineral_phase[i],name_of mineral_phase[j]); for (i1=1; il< =9; il++) fprintf (fo, "% 6.3f ", weight_of mineral_r[j][i1]); fprintf (fo, "\n"); for (i1=1; it < =9; il++) { fprintf (fo, " %6.3f ", weight_of mineral_r[i][il]); for (j1=1; jl<=9; j1++){ residual_square_array[i1][j1][I] = 0.0; for (k=1; k< =numberof oxide; k+ +){ sum = 0.0; for (k1=1; kl < =numberof mineral_phase; kl + +){ if (kl= =i) sum += di ff of m_and_d_save[k][kl]*weight_of mineral_r[kl][i1]; else { if (kl= =j) sum += di ff of m_and_d_save[k][kl]*weight_of mineral_r[kl] [j 1]; else sum + = diff of m_and_d_save[k][k1]*weight_of mineral_phase_r[kl]; } }  residual_r[k] = diff of_p_and_d[k] - sum; residual_square_array[i1][j1][I] + = residual_r[k]*residual_r[k]; }  residual_square_diff array[i1][j 1] [I] = residual_square_array[i1][j 1] [I]-residual_square; fprintf (fo, " %6.3f ", residual_square_diff array[il][jl][I]); 137  }  fprintf (fo, " \n"); }  }  }  /* 4./  /* calculation of the confidence level of the solutions around the best solution */ /* call gammq(a,x) function to calculate the delta_chisq on the given significant level p */ if (degrees_of freedom > 0){ i = 1; motification_of theleft_and_right_sides = 1.0; while (motification_of the_left_and_right_sides > 0){ left_side_of fitted_delta_chisq = i*0.001/2.0; left_side_of root = gammq(half degrees_of freedom, left_side_of fitted_delta_chisq) one_minus_p; right_side_of fitted_delta_chisq = (i + 1)*O. 001 /2.0; right_side_of root = gammq(half degrees_of freedom, right_side_of fitted_delta_chisq) one_minus_p; motification_of the_left_and_right_sides = left_side_of root * right_side_of root; i+ =1; }  printf ("motification_of the_left_and_right_sides = %20.18fin", motification_of theleft_and_right_sides); printf ("left_side_of fitted_delta_chisq = %f \n", left_side_of fitted_delta_chisq); printf ("right_side_of fitted_delta_chisq = %f \n", right_side_of fitted_delta_chisq); fitted_delta_chisq = (left_side_of fitted_delta_chisq + right_side_of fitted_delta_chisq); printf ("fitted_delta_chisq = %f \n", fitted_delta_chisq); alpha = sqrt(fitted_delta_chisq); printf ("alpha = square root of fitted delta chisq = %nn", alpha); fprintf (fo, " \n\n*****************************************************\ n ") ; fprintf (fo, "Del_Chisq Value = %f Less than Which the Probability of Chi-Sq Variable Is %f \n",^fitted_delta_chisq, p); fprintf (f0, "*****************************************************\ n\n\e) ; /* outputting the parameters of the solution confidence region */ fprintf (fo, " \n\nThe Orientations and Lengths of Semi-Axes of Chi-Squared Solution\n"); fprintf (fo," ^ \n\n"); for (i = 1; i< =number_of mineral_phase;i + +) {  fprintf (fo, "("); for (j = 1; j <number_of mineral_phase; j + +) fprintf (fo, " %5.3f, ", NUN); fprintf (fo, " %5.3f)^", v[number_of mineral_phase][i]); unit_length_of axis[i] = 1.0/w[i]; delta_chisq_length_of axis[i] = alpha*unit_length_of axis[i]; fprintf (fo, "Del-Chisq Length = %7.4f\n", delta_chisq_length_of axis[i]); 138  /* output the chi and residual square distribution for plotting profile into the data file plot.out */ /* output the chi_square distribution around the best solution */ if (numberof mineral_phase= =3){ fprintf (fo2, "\nChi_squared mapped against amounts of phases\n"); fprintf (fo2, " ^ \n"); I=0; for (i=1; i< =numberof mineral_phase; i+ +) for (j=i + 1; j < =number_of mineral_phase; j + +){ I+ =1; fprintf (fo2, " \n%s vs %s (", name_of mineral_phase[i], name_of_mineral_phase[j]); for (k=1; k< =numberof mineral_phase; k+ +){ if (k! =i) if (k!=j) fprintf (fo2, "%s = %7.4f gram ", name_of mineral_phase[k], weight_of mineral_phase[k]); }  fprintf (fo2, ")\n"); fprintf (fo2, " %s\\ %s ", name_of mineral_phase[i], name_of_mineral_phase[j]); for (i1=1; it < =9; il+ +) fprintf (fo2, " %6.3f ", weight_of mineral[j][il]); fprintf (fo2, "\n"); for (il =1; il < =9; il+ +) fprintf (fo2, " %6.3f^", weightof mineral[i][il]); for (j1=1; jl<=9; jl++) fprintf (fo2, " %6.2f ", chisq_array[il][jl][I]); fprintf (fo2, "\n"); {  }  }  /* output the distribution of the chi_square difference between the best solution and the chi_square around it */ fprintf (fo2, "\nDel-chi-squared against amounts of phases\n"); fprintf (fo2, " ^ \n"); I=0; for (i =1; i< =number of mineral_phase; i+ +) for (j=i +1; j < =numberof mineral_phase; j + +){ I+=1; fprintf (fo2, " \n%s vs %s (", nameof mineral_phase[i], name_of_mineral_phase[j]); for (k=1; k< =numberof mineral_phase; k+ +){ if (k!=i) if (k! =j) fprintf (fo2, " %s = %7.4f gram ", name_of mineral_phase[k], weight_of mineral_phase[k]); } fprintf (fo2, ")\n"); fprintf (fo2, " %s\\ %s ", nameof mineral_phase[i], name_of_mineral_phase[j]); for (i1=1; it < =9; il++) fprintf (fo2, " %6.3f ", weight_of mineral[j][il]); 139  fprintf (fo2, "\n"); for (i1=1; il < =9; il+ +) { fprintf (fo2, " %6.3f^", weight_of mineral[i][il]); for (j1=1; jl<=9; jl++) fprintf (fo2, " %6.2f ", chisq_diff array[il][j1][11); fprintf (fo2, "\n"); } } /* output the residual_square distribution around the best solution */ fprintf (fo2, "\nR-squared mapped against amounts of phases\n"); fprintf (fo2, " ^ \n"); I=0; for (i=1; i< =numberof mineral_phase; i+ +) for (j = i + 1; j < =numberof mineral_phase; j + +){ I+ =1; fprintf (fo2, " \n%s vs %s (", name_of mineral_phase[i], name_of mineral_phase[j]); for (k =1; k< =numberof mineral_phase; k+ +){ if (k!=i) if (k!=j) fprintf (fo2, "%s = %7.4f grams ", name_of mineral_phase[k], weight_of mineral_phase_r[k]);) fprintf (fo2, ")\n"); fprintf (fo2, " %s\\ %s ",name_of mineral_phase[i],name_of mineral_phase[j]); for (i1=1; it < =9; il+ +) fprintf (fo2, " %6.3f ", weight_of mineral_r[j][i1]); fprintf (fo2, "\n"); for (i1=1; it < =9; il+ +) { fprintf (fo2, " %6.3f ", weight_of mineral_r[i][il]); for (j1=1; jl<=9; jl++) fprintf (fo2, " % 6 . 3 f " , residual_square_array[il][jl][I]); fprintf (fo2, "\n"); }  }  /* output the distribution of the residual_square difference between the best solution and the residual_square around it */ fprintf (fo2, "\nDelta-R-squared mapped against amounts of phases\n"); fprintf (fo2, " ^ \n"); I=0; for (i =1; i < =number of mineral_phase; i+ +) for (j=i +1; j < =numberof mineral_phase; j + + )( I+ =1; fprintf (fo2, " \n %s vs %s (", name_of mineral_phase[i], name_of mineral_phase[j]); for (k =1; k< =numberof mineral_phase; k+ +){ if (k!=i) if (k!=j) fprintf (fo2, "%s = %7.4f grams ", name_of mineral_phase[k], weight_of mineral_phase_r[k]);} fprintf (fo2, ")\n"); 140  fprintf (fo2, " %s\ \ %s ",name_of mineral_phase[i],name_of mineral_phase[j]); for (i1=1; il< =9; il++) fprintf (fo2, " %6.3f ", weight_of mineral_r[j][il]); fprintf (fo2, "\n"); for (i1=1; il< =9; il++) fprintf (fo2, " %6.3f ", weight_of mineral_r[i][il]); for (j1=1; jl<=9;j1++) fprintf (fo2, " %6.3f ", residual_square_diff array[il][j1][I]); fprintf (fo2, "\n"); {  }  }  }  /* The SVDCMP.0 function called by main program. static float at,bt,ct; #define PYTHAG(a,b) ((at= fabs(a)) > (bt=fabs(b)) ? \ (ct=bt/at,at*sqrt(1.0+ct*ct)) : (bt ? (ct=at/bt,bt*sqrt(1.0 +ct*ct)): 0.0)) static float maxargl,maxarg2; #define MAX(a,b) (maxarg1=(a),maxarg2=(b),(maxargl) > (maxarg2) ?\ (maxargl) : (maxarg2)) #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a)) void svdcmp(a,m,n,w,v) float a[12][8], w[8], v[8][8]; int m,n; {  int flag,i,its,j,jj,k,l,nm; float c,f,h,s,x,y,z; float anorm=0.0,g=0.0,scale=0.0; float *rv1,*vector(); void nrerror(),free_vector(); if (m < n) nrerror("SVDCMP: You must augment A with extra zero rows"); rvl=vector(1,n); for (i=1;i< =n;i+ +) 1=i+1; rvl[i]=scale*g; g=s=scale=0.0; if (i < = m) for (k=i;k < =m;k+ +) scale + = fabs(a[k][i]); if (scale) { for (k=i;k<=m;k+ +) a[k][i] /= scale; s + = a[k][i]*a[k][i]; {  {  f=a[i][i]; g = -SIGN(sqrt(s),f); h=f*g-s;  141  a[i][i] = f-g; if (i != n) { for (j =1;j < =n;j + +) { for (s=0.0,k=i;k< =m;k+ +) s + = a[k][i]*a[k][j];  f=s/h; for (k=i;k< = m;k + +) a[k][j] + = f*a[k][i]; }  }  for (k=i;k< =m;k+ +) a[k][i] *= scale; }  w[i] = scale*g; g=s=scale=0.0; if (i < = m && != n) { for (k=1;k< =n;k+ +) scale + = fabs(a[i][k]); if (scale) { for (k=1;k<=n;k+ +) { a[i][k] /= scale; s + = a[i][k]*a[i][k]; }  f=a[i][1]; g = -SIGN(sqrt(s),0; h= f*g-s; a[i][1]=f-g; for (k=1;k< =n;k+ +) rvl[k]=a[i][k]/h; if (i != m) { for (j =1;j < =m;j+ +) { for (s=0.0,k=1;k< =n;k+ +) s += for (k=1;k< =n;k+ +) a[j][k] + = s*rvl[k]; }  }  for (k=1;k< =n;k+ +) a[i][k] *= scale; }  }  anorm=MAX(anorm,(fabs(w[i])+fabs(rvl[i]))); }  for (i=n;i > =1;i--) { if (i <n){ if (g) { for (j=1;j< =n;j+ +) v[j][i] = (a[i][j]/a[i][1])/g; for (j =1;j < =n;j+ +) { for (s =0.0,k =1;k < =n;k+ +) s + = a[i][k]*v[k][j]; for (k=1;k< =n;k+ +) v[k][j] + = s*v[k][i]; }  }  for (j----1;j < =n;j + +) v[i][j] =v[j][i] =0.0; -  }  v[i][i] =1.0; g=--rvl[i]; 1=--i; }  142  for (i=n;i > =1;i--) { 1=i+1; g=w[i]; if (i < n) for (j =1;j < =n;j+ +) a[i][j]=0.0; if (g) { g= 1.0/g; if (i != n) { for (j =1;j < =n;j+ +) { for (s=0.0,k=1;k< = m;k+ +) s += a[k][i]*a[k][j]; f=(s/a[i][i]rg; for (k=i;k< =m;k + +) a[k][j] + = f*a[k][i]; }  } for (j=i;j < =m;j+ +) a[j][i] *= g; } else { for (j=i;j < =m;j+ +) a[j][i]=0.0; } + +a[i][i]; } for (k=n;k > =1;k--) ( for (its= 1;its < =30;its+ +) { flag= 1; for (1=k;1> =1;1--) { nm=1-1; if (fabs(rv1[1])+anorm = = anorm) { flag=0; break; } if (fabs(w[nm])+anorm == anorm) break; } if (flag) { c=0.0; s= 1.0; for (i=1;i<=k;i+ +) { f=s*rvl[i]; if (fabs(f)+anorm != anorm) { g=w[i]; h=PYTHAG(f,g); w[i]=h; h= 1.0/h; c=g*h; s=(-f*h); for (j=1;j<=m;j+ +) { Y=alllinin]; z=a[j][i]; a[j][nm]=y*c+z*s; a[j][i]=z*c-y*s; }  }  } } z=w[k]; if (1 = = k) { 143  if (z < 0.0) { w[k] = -z; for (j=ld < = n;j + +) v[j][k]=(-v[j][k]); }  break; }  if (its = = 30) nrerror("No convergence in 30 SVDCMP iterations"); x=w[1]; nm=k-1; y =w[nm]; g=rvl[nm]; h=rvl[k]; f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y); g=PYTHAG(f,1.0); f=((x-z)*(x+z)+h*((y/(f+SIGN(g,f)))-h))/x; c=s=1.0; for (j =1;j < = mud + +) i=j+1; g=rvl[i]; Y= wEil; h=s*g; g=c*g; z=PYTHAG(f,h); rvl[j]=z; c= f/z; s=h/z; f=x*c+g*s; g=g*c-x*s; h= y*s; y=y*c; for (jj=1;jj < =ndj+ +) { x=v[jj][j]; z=v[jj][i]; v[jj][j]= x*c + z*s; v[j][i]=z*c-x*s; }  z=PYTHAG(f,h); w[j]=z; if (z) { z= 1.0/z; c=f*z; s=h*z; }  f=(c*g)+(s*y); x=(c*y)-(s*g); for (jj=1;jj < =m;jj+ +) {  Y=41.11U];  z= ii[jj][i]; a[jj][j]=y*c+z*s; a[jj][i]=z*c-y*s; }  rv1[1] =0.0; rv l[k] = f; 144  w[k]=x;  free_vector(rv1,1,n); }  #undef SIGN #undef MAX #undef PYTHAG  /* The SVBKSB.0 function called by main program. void svbksb(u,w,v,m,n,b,x) float u[12][8], w[8], v[8][8], b[S],x[8]; int m,n; {  int jj,j,i; float s,*tmp,*vector(); void free vector(); tmp = vector(1,n); for (j= 1;j < =n;j + +) { s =0.0; if (w{l]) { for (i = 1;i < =m;i+ +) s += u[i][j]*b[i]; s /= with }  tmp[j]=s; }  for (j =1;j < =n;j + +) { s=0.0; for (jj =1; jj < =n;jj+ +) s + = v[j][jj]*tmp[jj]; x[j]=s; }  free_vector(tmp,l,n); }  /* Function to count the number of characters in a string (using pointer). int string_length (string) char *string; {  char *cptr = string; while ( *cptr ) + +cptr; return (cptr - string); }  /* The Gammln.c function called by main program. ^  *1  float gammln(xx) float xx; {  double x,tmp,ser; static double cof[6]= {76.18009173,-86.50532033,24.01409822, -1.231739516,0.120858003e-2,-0.536382e-5}; int j; x =xx-1.0; tmp=x+ 5.5; tmp -= (x+0.5)*log(tmp); ser=1.0; for (j =0;j < =5;j+ +) { x + = 1.0; ser + = cof[j]/x; } return -tmp+log(2.50662827465*ser);  }  /* The Gammq.c function called by main program. ^  *1  float gammq(a,x) float a,x; { float gainser,gammcf,g1n; void gcf(),gser(),nrerror(); if (x < 0.0 I I a < = 0.0) nrerror("Invalid arguments in routine GAMMQ"); if (x < (a+1.0)) { gser(&gamser,a,x,&gln); return 1.0-gamser; } else { gcfi&gammcf,a,x,&gln); return gammcf; }  }  /* The Gcf.c function. ^  *1  #define ITMAX 100 #define EPS 3.0e-7 void gcfigammcf,a,x,g1n) float a,x,*gammcf,*gln; { int n; float gold=0.0,g,fac=1.0,b1=1.0; float b0=0.0,anf,ana,an,al,a0=1.0; float gammln(); void nrerror(); 146  *gln=gammln(a); al =x; for (n=1;n< =ITMAX;n+ +) { an= (float) n; ana = an-a; a0 = (al +a0*ana)*fac; b0 = (Ill +b0*ana)*fac; anf= an*fac; al =x*a0+anf*al; bl =x*b0+anf*bl; if (al) { fac =1.0/al ; g=bl*fac; if (fabs((g-gold)/g) < EPS) { *gammcf= exp(-x +a*log(x)-(*gln))*g; return; }  gold= g; }  }  nrerror("a too large, ITMAX too small in routine GCF"); }  #undef ITMAX #undef EPS /* The Gser.c function.  #define ITMAX 100 #define EPS 3.0e-7 void gser(gamser,a,x,g1n) float a,x,*gamser,*gln; {  int n; float sum,del,ap; float gammln(); void nrerror(); *gln=gammln(a); if (x < = 0.0) { if (x < 0.0) nrerror("x less than 0 in routine GSER"); *gamser = 0.0; return; } else { ap = a; del = sum= 1.0/a; for (n=1;n< =ITMAX;n+ +) { ap + = 1.0; del *= x/ap; sum + = del; if (fabs(del) < fabs(sum)*EPS) { 147  *gamser=sum*expex +a*log(x)-(*gln)); return; }  ) nrerror("a too large, ITMAX too small in routine GSER"); return; }  }  #undef ITMAX #undef EPS /* The Malloc.h function^ .  1* malloc.h: This file is referenced by the preprocessor directive #include <malloc.h> . You should either (i) copy this file to the appropriate directory so that it can be found, or (ii) replace the statements #include <malloc.h> by different #include statements, compatible with your compiler's handling of malloc() and free(). ^*1 char *malloc(); /* The Nrutil.c function .  void nrerror(error_text) char error_text[]; { void exit(); fprintfistderr,"Numerical Recipes run-time error... \n"); fprintfistderr," %s\n",error_text); fprintfistderr,"...now exiting to system... \n"); exit(1); }  float *vector(nl,nh) int nl,nh; { float *v; v= (float *)mallocaunsigned) (nh-n1+1)*sizeof(float)); if (!v) nrerror("allocation failure in vector()"); return v-nl; }  int *ivector(nl,nh) int nl,nh; { int *v;  v =(int *)malloc((unsigned) (nh-n1+1)*sizeof(int)); if (!v) nrerror("allocation failure in ivector()"); return v-nl; )  double *dvector(nl,nh) int nl,nh; {  double *v; v=(double *)mallocaunsigned) (nh-n1+1)*sizeof(double)); if (!v) nrerror("allocation failure in dvector()"); return v-nl;  )  float **matrix(nrl,nrh,ncl,nch) int nrl,nrh,ncl,nch; {  int i; float **m; m= (float **) malloc((unsigned) (nrh-nrl + 1)*sizeof(float*)); if (!m) nrerror("allocation failure 1 in matrix()"); m -= nrl; for(i=nr1;i < =nrh;i+ +) { m[i]=(float *) malloc((unsigned) (nch-nc1+1)*sizeof(float)); if (!m[i]) nrerror("allocation failure 2 in matrix()"); m[i] -= ncl; }  return m; }  double **dmatrix(nrl,nrh,ncl,nch) int nrl,nrh,ncl,nch; {  int i; double **m; m=(double **) malloc((unsigned) (nrh-nrl +1)*sizeof(double*)); if (!m) nrerror("allocation failure 1 in dmatrix()"); m -= nrl; for(i=nrl;i< =nrh;i+ +) { m[i] = (double *) malloc((unsigned) (nch-ncl + 1)*sizeof(double)); if (!m[i]) nrerror("allocation failure 2 in dmatrix()"); m[i] -= ncl; }  return m; }  int **imatrix(nrl,nrh,ncl,nch) 149  int nrl,nrh,ncl,nch; { int i,**m; m=(int **)malloc((unsigned) (nrh-nrI+1)*sizeof(int*)); if (!m) nrerror("allocation failure 1 in imatrix()"); m -= nrl; for(i=nrI;i< =nrh;i+ +) { m[i]= (int *)malloc((unsigned) (nch-ncl + 1)*sizeof(int)); if (!m[i]) nrerror("allocation failure 2 in imatrix0"); m[i] -= ncl; } return m; }  float **submatrix(a,oldrl,oldrh,oldcl,oldch,newrl,newc1) float **a; int oldrl,oldrh,oldcl,oldch,newrl,newcl; {  int i,j; float **m; m= (float **) malloc((unsigned) (oldrh-oldr1+1)*sizeof(float*)); if (!m) nrerror("allocation failure in submatrix()"); m -= newrl; for(i=oldrl,j =newr1;i < =oldrh;i+ +,j + +) m[j]=a[i] +oldcl-newcl; return m;  }  void free_vector(v,n1,nh) float *v; int nl,nh; {  free((char*) (v +n1));  }  void free_ivector(v,n1,nh) int *v,n1,nh; { free((char*) (v +n1)); }  void free_dvector(v,n1,nh) double *v; int nl,nh; { free((char*) (v +n1)); 150  }  void free_matrix(m,nrl,nrh,ncl,nch) float **m; int nrl,nrh,ncl,nch; {  int i; for(i=nrh;i > =nrl;i--) free((char*) (m[i]+ncl)); free((char*) (m+nr1)); }  void free_ dmatrix(m,nrl,nrh,ncl,nch) double **m; int nrl,nrh,ncl,nch; int i;  }  for(i=nrh;i > =nrl;i--) free((char*) (m[i]+ncl)); free((char*) (m+nr1));  void free_imatrix(m,nrl,nrh,ncl,nch) int **m; int nrl,nrh,ncl,nch;  {  }  int i; for(i=nrh;i > =nrl;i--) free((char*) (m[i]+ncl)); free((char*) (m+nr1));  void freesubmatrix(b,nrl,nrh,ncl,nch) float **b; int nrl,nrh,ncl,nch;  {  }  free((char*) (b+nr1));  float **convert_matrix(a,nrl,nrh,ncl,nch) float *a; int nrl,nrh,ncl,nch;  {  int i,j,nrow,ncol; float **m; nrow=nrh-nr1 +1; ncol = nch-ncl +1; 151  m = (float **) malloc((unsigned) (nrow)*sizeof(float*)); if (!m) nrerror("allocation failure in convert_matrix0"); m -= nrl; for(i =0,j =nr1;i < = nrow-1 ;i + + ,j + + ) m[j] =a +ncol*i-ncl; return m; }  void freeconvert_matrix(b,nrl,nrh,ncl ,nch) float **b; int nrl,nrh,ncl,nch; { }  free((char*) (b + nrl));  /* The Nrutil.h function.  float *vector(); float **matrix(); float **convert_matrix(); double *dvector(); double **dmatrix(); int *ivector(); int **imatrix(); float **submatrix(); void free_vector(); void free_dvector(); void free_ivector(); void free_matrix(); void free_dmatrix(); void freeimatrix(); void free_submatrix(); void free_convert_matrix(); void nrerror();  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0052616/manifest

Comment

Related Items