On the misuse of regression in earth science

Simple least-squares regression analysis is applied to almost all empirical curve-fitting problems in earth science (and related fields). Its use, however, should be restricted to predictive situations. For comparisons with theory or among fitted lines, the related technique termed functional analysisshould be employed. To apply this method, the ratio of the random components of the variances of the variables must be estimated. Principles are illustrated with examples from geomorphometry, especially the stream frequency-drainage density relation.     

A pseudo-genetic algorithm suited to compositional data

Genetic algorithms can solve least-squares problems where local minima may trap more traditional methods. Although genetic algorithms are applicable to compositional as well as noncompositional data, the standard implementation treats compositional data awkwardly. A need to decode, renormalize, then reincode the fitted parameters to regain a composition is not only computationally costly, but may thwart convergence. A modification to the genetic algorithm, described here, adapts the tools of reproduction, crossover, and mutation to compositional data. The modification consists of replacing crossover with a linear mixture of two parents and replacing mutation with a linear mixture of one of the members of the breeding population and a randomly generated individual. By using continuously evolving populations, rather than discrete generations, reproduction is no longer required. As a test of this new approach, a mixture of four Gaussian functions with given means and variances are deconvolved to recover their mixing proportions.     

The New “York” Regression: Application of an Improved Statistical Method to Geochemistry

York's (1969) method of regression, determining the best-fit line to data with errors in both variables using a least-squares solution, has become an integral part of isotope geochemistry. Although other methods agree with York's best-fit line (e.g., maximum likelihood), there is little agreement on the standard-error estimates for slope and intercept values. The reasons for this are differing levels of approximation used to compute the standard error, doubts concerning procedures for determining a confidence interval once the standard error has been estimated, and a typographical error in the original publication. This paper examines York's method of regression and standard errors of the parameters of a best-fit line. A very accurate method for determining the standard error in slope and intercept values is introduced, which eliminates the need to multiply the standard-error estimate by the goodness-of-fit parameter known as MSWD. In addition, a derivation of a fixed-intercept method of regression is introduced, and interpretations of MSWD and use of the t-adjustment in confidence intervals are discussed. The accuracy of the standard-error computations is determined by comparing the results to slope and intercept statistics generated from several thousand Monte Carlo regressions using synthetic ⁴⁰Ar/³⁹Ar inverse isochron data.     

An Empirical Method for Optimal Robust Regional-Residual Separation of Geophysical Data

Prior to interpretation and further analysis, many datasets must first be separated into regional and residual components. Traditional techniques are either subjective (e.g., graphical methods) or nonrobust (e.g., all least-squares based methods). Bathymetric data, with their broad spectrum, pose serious difficulties to these traditional methods, in particular those based on spectral decomposition. Spatial median filters offer a solution that is robust, objective, and often defines regional components similar to those produced graphically by hand. Characteristics of spatial median filters in general are discussed and a new empirical method is presented for determining the width of the robust median filter that accomplishes an optimal separation of a gridded dataset into its regional and residual components. The method involves tracing the zero-contour of the residual component and evaluating the ratio between the volume enclosed by the surface inside this contour and the contour's area. The filter width giving the highest ratio (or mean amplitude) is called the Optimal Robust Separator (ORS) and is selected as the optimal filter producing the best separation. The technique allows a unique and objective determination of the regional field and enables researchers to undertake reproducible separations of regional and residual components. The ORS method is applied to both synthetic data and bathymetry/topography of the Hawaiian Islands; ways to improve the technique using alternative diagnostic quantities are discussed.     

Performance studies on six heuristic global optimization methods in the location of critical slip surface

Recently, many heuristic global optimization algorithms have evolved with success for treating various types of problems. Majority of these algorithms have not been applied to slope stability problem for which the presence of soft band and convergence problem (discontinuity of the objective function) may create difficulties in the minimization process. In this paper, six heuristic optimization algorithms are applied to some simple and complicated slopes. The effectiveness and efficiency of these algorithms under different cases are evaluated, and it is found that no single method can outperform all the other methods under all cases, as different method has different behavior in different types of problems. For normal cases, the particle swarm method appears to be effective and efficient over various conditions, and this method is recommended to be used. For special cases where the objective function is highly discontinuous, the simulated annealing method appears to be a more stable solution.     

Nonlinear solver for three-phase transport problems based on approximate trust regions

Implicit transport solvers used in reservoir simulation can take longer time steps than explicit solvers, but for long time steps, the commonly used Newton-Raphson’s method will often fail to converge. The convergence issues may manifest themselves as oscillating residuals even though the implicit discretization itself is stable. This behavior occurs because the fractional flow-type flux functions often change between convex and concave during long time steps, resulting in multiple contraction regions for the Newton-Raphson solver. The common strategy to overcome this is to set limits on the saturation changes during the nonlinear iteration, but such a limit has to be determined on a case by case basis, excess iterations may be required, and practical convergence is not guaranteed for a given problem. Previous work on this problem by multiple authors has resulted in solvers based on trust regions, where unconditional convergence can be obtained for incompressible two-phase flow provided a priori analytical knowledge of the flux function exists. The goal of our work is to extend this methodology to a solver where inflection points demarking the different contraction regions do not need to be explicitly known. Instead, these values are estimated during the solution process, giving improved convergence by a local computation for each interface in the simulation model. By systematically reducing updates over regions known to produce convergence issues, it is possible to greatly reduce the computational expense, making the same formulation suitable for an arbitrary number of components. We present a series of numerical results, including arbitrary time-step lengths for two and three-phase gravity segregation, as well as three-dimensional gas and water injection problems with wells and a mixture of both viscous and gravity-dominated flow regimes. The test cases are a systematic validation on a wide variety of both analytical and tabulated relative permeability curves.     

Numerical manifold method for dynamic nonlinear analysis of saturated porous media

It is well known that the Babuska–Brezzi stability criterion or the Zienkiewicz–Taylor patch test precludes the use of the finite elements with the same low order of interpolation for displacement and pore pressure in the nearly incompressible and undrained cases, unless some stabilization techniques are introduced for dynamic analysis of saturated porous medium where coupling occurs between the displacement of solid skeleton and pore pressure. The numerical manifold method (NMM), where the interpolation of displacement and pressure can be determined independently in an element for the solution of u–p formulation, is derived based on triangular mesh for the requirement of high accurate calculations from practical applications in the dynamic analysis of saturated porous materials. The matrices of equilibrium equations for the second‐order displacement and the first‐order pressure manifold method are given in detail for program coding. By close comparison with widely used finite element method, the NMM presents good stability for the coupling problems, particularly in the nearly incompressible and undrained cases. Numerical examples are given to illustrate the validity and stability of the manifold element developed. Copyright © 2006 John Wiley & Sons, Ltd.     

Fast automatic solution of the inverse resistivity problem

The paper presents a fast automatic approach to solve the inverse resistivity problem, assisted by optimization, which is a non-linear model-fitting technique. The selected inverse problems are ill-posed and the inverse solution is defined by ‘best fit’ in the sense of least-squares. Formulations are presented in a systematic manner for Newton’s method, least squares method and Marquardt’s modification (ridge regression) method based on local linearization of non-linear problem. The convergence of least-squares method and Marquardt’s method, to provide a robust solution, are first tested on a theoretical model and effectiveness of Marquardt’s method is demonstrated, and then two-field apparent resistivity curves from Banda district, India are interpreted and discussed.     

Error propagation in equations for geochemical modeling of radiogenic isotopes in two-component mixing

This paper presents error propagation equations for modeling of radiogenic isotopes during mixing of two components or end-members. These equations can be used to estimate errors on an isotopic ratio in the mixture of two components, as a function of the analytical errors or the total errors of geological field sampling and analytical errors. Two typical cases (“Small errors” and “Large errors”) are illustrated for mixing of Sr isotopes. Similar examples can be formulated for the other radiogenic isotopic ratios. Actual isotopic data for sediment and basalt samples from the Cocos plate are also included to further illustrate the use of these equations. The isotopic compositions of the predicted mixtures can be used to constrain the origin of magmas in the central part of the Mexican Volcanic Belt. These examples show the need of high quality experimental data for them to be useful in geochemical modeling of magmatic processes.     

Multiple component end-member mixing model of dilution: hydrochemical effects of construction water at Yucca Mountain, Nevada, USA

The standard dual-component and two-member linear mixing model is often used to quantify water mixing of different sources. However, it is no longer applicable whenever actual mixture concentrations are not exactly known because of dilution. For example, low-water-content (low-porosity) rock samples are leached for pore-water chemical compositions, which therefore are diluted in the leachates. A multicomponent, two-member mixing model of dilution has been developed to quantify mixing of water sources and multiple chemical components experiencing dilution in leaching. This extended mixing model was used to quantify fracture-matrix interaction in construction-water migration tests along the Exploratory Studies Facility (ESF) tunnel at Yucca Mountain, Nevada, USA. The model effectively recovers the spatial distribution of water and chemical compositions released from the construction water, and provides invaluable data on the matrix fracture interaction. The methodology and formulations described here are applicable to many sorts of mixing-dilution problems, including dilution in petroleum reservoirs, hydrospheres, chemical constituents in rocks and minerals, monitoring of drilling fluids, and leaching, as well as to environmental science studies.     

Thermochemical parameters of minerals from oxygen-buffered hydrothermal equilibrium data: Method,application to annite and almandine

Reversed univariant hydrothermal phase-equilibrium reactions, in which a redox reaction occurs and is controlled by oxygen buffers, can be used to extract thermochemical data on minerals. The dominant gaseous species present, even for relatively oxidizing buffers such as the QFM buffer, are H₂O and H₂; the main problem is to calculate the chemical potentials of these components in a binary mixture. The mixing of these two species in the gas phase was assumed by Eugster and Wones (1962) to be ideal; this assumption allows calculation of the chemical potentials of the two components in a binary gas mixture, using data in the literature. A simple-mixture model of nonideal mixing, such as that proposed by Shaw (1967), can also be combined with the equations of state for oxygen buffers to permit derivation of the chemical potentials of the two components. The two mixing models yield closely comparable results for the more oxidizing buffers such as the QFM buffer. For reducing buffers such as IQF, the nonideal-mixing correction can be significant and the Shaw model is better.The procedure of calculation of mineralogical thermochemical data, in reactions where hydrogen and H₂O simultaneously appear, is applied to the experimental data on annite, given by Wones et al. (1971), and on almandine, given by Hsu (1968). For annite the results are: Standard entropy of formation from the elements, S _f ⁰ (298, 1)=–283.35±2.2 gb/gf, S ⁰ (298, 1) =+92.5 gb/gf. G _f ⁰ (298, 1)=–1148.2±6 kcal, and H _f ⁰ (298, 1)=–1232.7±7 kcal. For almandine, the calculation takes into account the mutual solution of FeAl₂O₄ (Hc) in magnetite and of Fe₃O₄ (Mt) in hercynite and the temperature dependence of this solid solution, as given by Turnock and Eugster (1962); the calculations assume a regular-solution model for this binary spinel system. The standard entropy of formation of almandine, S _f,A ⁰ (298, 1) is –272.33±3 gb/gf. The third law entropy, S ⁰ (298, 1) is +68.3±3 gb/gf, a value much less than the oxide-sum estimate but the deviation is nearly the same as that of grossularite, referring to a comparable set of oxide standard states. The Gibbs free energy G _f,A ⁰ (298, 1) is –1192.36±4 kcal, and the enthalpy H _f,A ⁰ (298, 1) is –1273.56±5 kcal.Publication authorized by the Director, U. S. Geological Survey.     

Robust statistics and geochemical data analysis

Advantages of robust procedures over ordinary least-squares procedures in geochemical data analysis is demonstrated using NURE data from the Hot Springs Quadrangle, South Dakota, U.S.A. Robust principal components analysis with 5% multivariate trimming successfully guarded the analysis against perturbations by outliers and increased the number of interpretable factors. Regression with SINE estimates significantly increased the goodness-of-fit of the regression and improved the correspondence of delineated anomalies with known uranium prospects. Because of the ubiquitous existence of outliers in geochemical data, robust statistical procedures are suggested as routine procedures to replace ordinary least-squares procedures.     

沉积岩系互稀释混合元素丰度关系

二元混合成国岩石元素丰度（ｘ，ｙ）遵守分式线性函数定律:式中ａ，ｂ，ｃ，ｄ为丰度关系常数，它们只和岩系中两个已知岩石（１和２）中元素丰度（ｘ_１，ｙ_１；ｘ_２，ｙ_２）有关，并等于其镜像∑变换（混合变换）（ａ，ｂ，ｃ，ｄ）＝∑∑|ｘ_１，ｙ_１；ｘ_２，ｙ_２|（２）元素丰度关系常数中，如果ｂ与Ｃ之比值小于零时，元素丰度ｙ＝Ｘ呈互稀释混合关系。     

The energetics of natural garnet solid solution

Approximate mixing properties of the end-member components of the quarternary garnet solid solution, (Fe,Mg,Ca,Mn)₃Al₂Si₃O₁₂, have been derived through theoretical analysis of observational data, combined with certain experimental results and crystal chemical considerations. The results suggest that the mixing of pyrope with grossularite, spessartite, and almandine would involve significant positive excess free energies of mixing leading to the critical mixing temperatures of 694±55, 535±140, and 479±63 °C respectively. Spessartite would mix with almandine nearly ideally, and with grossularite with small positive deviation from ideality. The quarternary solution reduces essentially to a ternary mixture of pyrope, grossularite, and almandine + spessartite. The solid solubility relation, and tie line coordinates in this ternary system has been calculated as a function of temperature; the solid solution is found to be intrinsically stable for practically all ternary compositions at 600 °C.     

Difference between model parameter estimation and curve fitting

Within a single well, a least-squares model such as Ø = a +bT _sonic may be appropriate, insofar as random errors in core measurement and in sonic tool response are concerned. However, if several wells are in the reservoir, and if structural variation across the field exists, no single least-squares solution will adequately describe all wells with the same precision. Interpolation, via a cubic spline, is offered as an alternative to model this type of variation. However, to control both inter- and intra-well variability, a combination of least squares and interpolation is the model of choice.     

Degree of Geochemical Similarity (DOGS): a simple statistical method to quantify and map affinity between samples from multi-element geochemical data sets

The 42 element, 1190 sample Mobile Metal Ion subset of the National Geochemical Survey of Australia database was used to develop and illustrate the concept of Degree of Geochemical Similarity of soil samples. Element concentrations were unified to parts per million units and log(10)-transformed. The degree of similarity of pairs of samples of known provenance in the Yilgarn Craton was obtained using least-squares linear regression analysis and demonstrated that the method successfully assessed the degree of similarity of soils related to granitoid and greenstone lithologies. Exploratory Data Analysis symbol maps of all remaining samples in the database against various reference samples were used to obtain correlation maps not only for granitoid- and greenstone-related soil types but also to distinguish between, for example, samples derived from marine vs regolith carbonate. Likewise, the distribution of soil samples having a geochemical fingerprint similar to mineralised provinces (e.g. Mount Isa) can be mapped, and this can be used as a first-order prospection tool.     

Solving the Buckley–Leverett equation with gravity in a heterogeneous porous medium

Immiscible two‐phase flow in porous media can be described by the fractional flow model. If capillary forces are neglected, then the saturation equation is a non‐linear hyperbolic conservation law, known as the Buckley–Leverett equation. This equation can be numerically solved by the method of Godunov, in which the saturation is computed from the solution of Riemann problems at cell interfaces. At a discontinuity of permeability this solution has to be constructed from two flux functions. In order to determine a unique solution an entropy inequality is needed. In this article an entropy inequality is derived from a regularisation procedure, where the physical capillary pressure term is added to the Buckley‐Leverett equation. This entropy inequality determines unique solutions of Riemann problems for all initial conditions. It leads to a simple recipe for the computation of interface fluxes for the method of Godunov. This revised version was published online in July 2006 with corrections to the Cover Date.     

End-member modeling of compositional data: Numerical-statistical algorithms for solving the explicit mixing problem

Linear mixing models of compositional data have been developed in various branches of the earth sciences (e.g., geochemistry, petrology, mineralogy, sedimentology) for the purpose of summarizing variation among a series of observations in terms of proportional contributions of (theoretical) end members. Methods of parameter estimation range from relatively straightforward normative partitioning by (nonnegative) least squares, to more sophisticated bilinear inversion techniques. Solving the bilinear mixing problem involves the estimation of both mixing proportions and end-member compositions from the data. Normative partitioning, also known as linear unmixing, thus can be regarded as a special situation of bilinear unmixing with (supposedly) known end members. Previous attempts to model linear mixing processes are reviewed briefly, and a new iterative strategy for solving the bilinear problem is developed. This end-member modeling algorithm is more robust and has better convergence properties than previously proposed numerical schemes. The bilinear unmixing solution is intrinsically nonunique, unless additional constraints on the model parameters are introduced. In situations where no a priori knowledge is available, the concept of an “ optimal ” solution may be used. This concept is based on the trade-off between mathematical and geological feasibility, two seemingly contradictory but equally desirable requirements of the unmixing solution.     

An interactive global optimization algorithm for geological problems

A number of problems in geology can be formulated so that they consist of optimizing a real-valued function (termed the objective function) on some interval or over some region. Many methods are available for solution if the function is unimodal within the domain of interest. Direct methods, involving only function evaluations, are particularly useful in geological problems where the objective function may be strongly nonlinear and constructed from sampled data. In practical problems, the objective function often is not unimodal. Standard optimization routines are not capable of distinguishing between local extrema or of locating the global extremum, which is the point of interest in most cases. The usual approach—trying several different starting points in the hope that the best local extremum found is the global extremum—is inefficient and unreliable. An ancillary algorithm has been developed which avoids these problems and which couples with a variety of local optimization routines. The algorithm first constructs a grid of objective function values over some feasible region. The region dimensions and grid spacings are based on specific problem considerations. First differences are then calculated for successive points along each grid line and monitored in sign only, which rapidly locates extrema. User interaction determines how many of these extrema will undergo further investigation, which is carried out by passing locations to a local optimization subroutine. The algorithm has proved successful on a number of problems. A geological example—determination of benthic mixing parameters in deep-sea sediments via minimization of stratigraphic offset between ¹⁸ O signals from two different species of planktonic foraminifera—is given. FORTRAN code is provided for the global optimization routine, a golden section search subroutine for one-dimensional objective functions, and a simplex subroutine for multidimensional problems.