New Techniques on Tensor Decomposition and Multidimensional Forward Modeling Chairpersons: J. Ledo and U. Schmucker 4.1 ALL-TIME APPARENT RESISTIVITY FOR TRANSIENT ELECTROMAGNETIC METHOD Bai Denghai(1) and Maxwell A. Meju(2) (1) Institute of Geology, China Seismological Bureau, Beijing 100029, P. R. China (2) Department of Geology, University of Leicester, Leicester LE1 7RH, United Kingdom bdh607@sina.com Apparent resistivity is an important response parameter in the analysis of transient electromagnetic soundings. In routine studies, this is commonly obtained using the early-time or late-time approximations leading to significant errors in the intermediate time range. The all-time apparent resistivity was suggested as a solution to this problem and there are two main approaches to the computation. One technique uses the combined early-time and late-time approximations to generate a single interpretive sounding curve but the apparent resistivity for intermediate time is poorly defined in this case. Another approach directly calculates apparent resistivity by application of exact formulae. Routine implementation of these schemes require significant numerical skill. Computation of all-time apparent resistivities based on the iterative method and the series expansion method may sometimes yield meaningless results for the intermediate time range for some particular cases, and has even led some workers to suggest that the early-time apparent resisitivities might be affected by the so-called "static shift''. A new method for calculating all-time apparent resistivity sounding curves for time-domain transient electromagnetic methods is proposed in this paper. Using this method, the all-time apparent resistivities can be uniquely defined. We found that the resulting apparent resistivities are more accurate than those based on early-time and late-time definitions for the central loop configuration. We have not found any evidence of static shift using this approach. 4.2 COMPARATIVE STUDY OF 3D FINITE-DIFFERENCE AND INTEGRAL EQUATION METHODS (COMMEMI-3 MODEL SET) Ivan Varentsov(1), Elena Fomenko(1), Nikolay Golubev(2), Salah Mehanee(2), Gabor Hursan(2), and Michael Zhdanov(2) (1) Geoelectromagnetic Research Institute RAS, 142090, P.O.B. 30, Troitsk, Moscow Region, Russia (2) University of Utah, Salt Lake City, UT 84112, USA igemi1@pop.transit.ru We held the comparative study of 3D finite-difference (FD) and integral equation (IE) codes, available at the University of Utah and at the Geoelectromagnetic Research Institute. A new set of test models, COMMEMI-3, was considered, containing two old COMMEMI 3-D models (Zhdanov et al., 1997), 3D1 and 3D2; and a new model, 3D3, combining several anomalies of different sizes and depths with high conductivity contrast (up to 30000). The simplest model, 3D1, was extended to include two variants of conductivity contrast (200 and 10000). Another high contrast model, suggested by Mackie (1999), was also studied. The IE methods are represented by SYSEM solution (Xiong, 1992) and an approximate quasi-analytical (QA) series (Zhdanov, Fang, 1997; Zhdanov et al, 1999). The common basis for two considered FD techniques is the application of the staggered grids, the use of modern conjugate gradient and Chebyshev iterative solvers with proper preconditioning, and monitoring of the EM field divergence for accuracy control. The first technique implements the primary total magnetic field FD solution (Varentsov, 1999), while the second deals with the anomalous electric field FD problem (Fomenko, 1999). The joint use of both techniques gives an independent error control at the FD approximation stage, and the fast, stable and almost monotonously convergent FD solvers give a new quality at the stage of the FD system solution. All modeling codes demonstrated sufficient accuracy and computational effectiveness for all examined models. The FD techniques have shown high accuracy and stability of results, quite comparable with the traditional accuracy level of IE methods. Both of them provided close results in spite of the difference in the formulation of initial boundary value problems. We observed a good convergence of the QA series in complicated applications. The SYSEM method is not as fast as other techniques in the case, when the discrete system is too large to be solved be the direct solver. The insufficient accuracy of the SYSEM results at some vertical conductivity interfaces are related to discretization problems and slow convergence of block-iteration scheme. 4.3 AN APPLICATION OF TENSOR INVARIANT ANALYSIS TO REAL DATA A.K. Agarwal(1), J.T. Weaver(1), and F.E.M. Lilley(2) (1) School of Earth and Ocean Sciences, University of Victoria, Victoria, BC, Canada (2) Research School of Earth Sciences, The Australian National University, Canberra, ACT, Australia numod@uvvm.uvic.ca weaver@phys.uvic.ca At the 14th Workshop we introduced seven independent rotational invariants of the magnetotelluric impedance tensor, plus one additional dependent invariant, all of which had a certain physical interpretation when becoming vanishingly small. The seventh independent invariant and the one dependent invariant have since been slightly modified for greater clarity in interpretation, but without affecting their general properties. In this paper we analyse three MT datasets from Australia, Germany and Canada respectively, by investigating the behaviour of the seven invariants of the measured impedances at the sites along each profile of recording stations. Colour-coded diagrams, classifying the various physical interpretations that can arise, are constructed with the sites positioned horizontally and with increasing period displayed vertically, so that regions of one-, two- and three-dimensionality are immediately apparent. In those regions where the structure is determined to be regionally two-dimensional--with or without small- scale distortion--the tensor invariants are used to recover a strike angle as well as the principal impedances. Results will be presented for all three datasets. 4.4 SCATTERING ON JOINT ELECTRICAL AND MAGNETIC INHOMOGENEITIES: OVERVIEW OF NONLINEAR APPROXIMATIONS, COMPARISON, LOGGING APPLICATIONS Arvidas Cheryauka and Michael Zhdanov Consortium for EM modeling and inversion (CEMI), University of Utah, USA acher@mines.utah.edu We extend linear and nonlinear approximations for electromagnetic fields in a medium with inhomogeneous distribution of both electrical and magnetic material properties. These approximations are presented in the form of tensor integrals over a domain with anomalous parameters. The developed approximations combine the linear and nonlinear estimations depending on the ratio of complex conductivity and magnetic susceptibility perturbations. These approximations form a basis for fast EM modeling and imaging in multi-dimensional environments where joint electrical and magnetic inhomogeneity is an essential feature of the model. Numerical tests carried out for one-dimensional electromagnetic logging applications demonstrate the validity of the theory and the effectiveness of the proposed approximations. 4.5 COMPARATIVE STUDY OF RELAXATION METHODS FOR EM INDUCTION IN HETEROGENEOUS SPHERES: MRA, BICG AND BICGSTAB Hiroaki Toh(1), Adam Schultz(2), and Makoto Uyeshima(3) (1) Dept Earth Sciences, Toyama University, Toyama 9308555, Japan (2) Earth Sciences Department, University of Wales, Cardiff CF1 3YE UK (3) Earthquake Research Institute, University of Tokyo, Tokyo 1130032, Japan toh@gdyn01.sci.toyama-u.ac.jp EM induction in fully heterogeneous spheres was solved by a finite difference method on a staggered grid. The forward problem thus formulated results in inverting a large sparse complex symmetric matrix. Since the matrix inversion can now be performed by iterative relaxation methods such as conjugate gradient methods, the problem was first solved by Mackie and Madden's (1993) minimum residual accelerated (MRA) method. However, since it is preferable to make the coefficient matrix Hermitian for robust application of the forward solver to any form of three-dimensional inversions, we rewrote the solver into a biconjugate gradient (BiCG) method valid even for complex symmetric matrices. A stabilized variant of the BiCG method, BiCGSTAB, was further employed to reduce the number of iterations per relaxation. Comparison of the three solvers showed that they had almost the same accuracy while the BiCGSTAB method was the fastest. The BiCGSTAB solver was also optimized over multi-frequencies. It was shown that substitution of the solution of the previous frequency as an intial solution for the subsequent frequency greatly improved the computation speed. It was also found that to sort the frequencies in ascending order yielded better speed and accuracy. In addition, the BiCGSTAB solver was further accelerated by frequency-by-frequency parallelization. Three-dimensional inversions of the global electrical conductivity in the mantle using the BiCGSTAB solver are now quite feasible since the combination of the BiCGSTAB relaxation, the initialization by the previous solution, the ascending sorting and the frequency-by-frequency parallelization gives an order of magnitude acceleration in total. A very early results of spherical inversions using this solver will be included in the presentation. 4.6 EXTENDED DECOMPOSITION OF MT DATA Xavier Garcia and Alan G. Jones Geological Survey of Canada, 615 Booth St., Ottawa, Ontario, K1A 0E9, Canada garcia@cg.NRCan.gc.ca Decomposition of magnetotelluric data into local galvanic 3D distortion matrix and a regional 2D Earth has caused a quantum leap in our understanding of complex data and our ability to handle those data. The Groom-Bailey method is the most widely adopted tensor decomposition approach, and rightly so given its physical basis and its separation of distortion parameters into determinable and indeterminable parts. However, on occasion this 3D over 2D, 3D/2D, decomposition fails in that the misfit of the model to the data is far greater than the data errors permit, and this failure is due either to inappropriately small error estimates for the data, or to the model being invalid. In this work we describe and demonstrate our attempts to extend MT tensor decomposition to 3D distortion of regional 3D data, 3D/3D. There are insufficient data to accomplish this uniquely for a single MT site, so some approximations must be made. The approach we use is to assume that two neighboring sites sense the same regional structure if they are sufficiently close compared to the skin depth to the structure, but that the two sites have differing galvanic distortion matrices. 4.7p CONTROLLED SOURCE MAGNETOTELLURICS: THE ROLE OF THE MAGNETIC FIELD H.M. Bibby(1), T.G. Caldwell(2) (1) Institute of Geological and Nuclear Sciences, PO Box 30368, Lower Hutt, New Zealand (2) Currently at Applied Geophysics Unit, University College Galway, Galway, Ireland h.bibby@gns.cri.nz Previous work has shown that there are considerable advantages to be gained from the analysis of long-offset time-domain electromagnetic data using the instantaneous apparent resistivity tensor. The tensor is defined by the relationship between the measured electric field at any instant, t, and the (time independent) current density in a uniform half space, ie. E(t) = Rho(t) J . In this paper we explore the equivalent relationship in the frequency domain, E(w) = Rho(w) J , and its application in Controlled Source Magnetotellurics (CSMT). This relationship can be used to define a method by which the magnetotelluric impedance tensor can be separated into electric and magnetic components each of which can be used to independently define a tensor apparent resistivity. The frequency domain apparent resistivity tensor Rho(w) is dependent on the electric field alone and is defined for all values of w. The limit reached by Rho(w) as w tends to infinity is the same as that of the equivalent limit of Rho(t) in the time domain (t tends to 0). Similarly at the other extreme as w tends to 0, Rho(w) tends to the dc apparent resistivity tensor (equal to Rho(t) as t tends to infinity). The apparent resistivity derived from the magnetic field alone is valid only in the far field zone (r/delta) >> 1, where the MT approximation for CSMT is applicable. This approach provides a number of insights into the characteristics of both MT and CSMT. 4.8p AMT DATA PROCESSING: SOME DIFFERENCE COMPARED TO MT ONE Anton Kocherov St. Petersburg State University, Universitetskaya emb. 7/9, St. Petersburg, 199034, Russia elmag@elmag.ecri.pu.ru The EM field spectrum observable in AMT frequency range reveals very often a presence of numerous sinusoidal components of high intensity at the frequency of power line and its higher harmonics up to frequencies of several kilohertz, when making surveys in highly industrialised areas. Unlike the MT-case, when peaks in the spectrum arise from geomagnetic pulsations used in processing together with noise-like intervals of records, in the AMT-case such peaks are related to industrial noise and the estimates of impedance at these frequencies are strongly biased from the "right" sounding curve obtained using the natural field. To be able to exclude such data, a spectral resolution must be sufficient to resolve peaks of adjoining harmonics and obtain some points in the spectrum between them related to the undistorted natural field. Two main parameters control the resolution in spectral estimation. The first is the number of data points used in FFT. Though it varies in MT usually from 32 to 1024, much more points may be necessary in AMT processing. The second parameter of interest is the shape of the "tapering" window used before FFT. As industrial fields can be several orders greater than the natural field in the magnitude, the side lobes aside such an intensive peak can completely distort spectral estimates at the "natural field" frequencies all over between the adjoining peaks. A Nuttall's window having extremely small side lobes is recommended. Frequently used averaging impedance values over groups of adjoining frequencies can be implied only after excluding distorted data. Finally, most of intensive atmospherics would be rejected in some processing techniques because of their dominating linear polarisation, though they show a higher signal to noise ratio and yield valid impedance estimates. 4.9p TENSOR DECOMPOSITION EXAMPLES BY USING THE CORRECTED PHASE DEVIATION FORMULAS Erno Pracser(1) and Laszlo Szarka (1) Lorand Eotvos Geophysical Institute, H-1145 Budapest, Kolumbusz u. 17-23, Hungary (2) Geodetic and Geophysical Researh Institute of the Hungarian Academy of Sciences, H-9401 Sopron POB 5, Hungary pracser@elgi.hu Electromagnetic field distortion due to elongated near-surface bodies is almost independent of the direction of the inducing field as it was numerically studied in details by Szarka L., Menvielle M., Tarits P., Adam A. 1994, Acta Geod. Geoph. Hung., 29, 125-138. The disturbing field depends from the position, orientation and conductivity of the disturbing body. Consequently, the distortion field in such cases is easily obtained from the static field due to an electric dipole. In a recent paper by Pracser E., Szarka L. 1999, Earth Planets Space, 51, 1019-1022 it has been shown that a correction of the phase deviation method for tensor decomposition provides more understandable and interpretable results than the original formulas published in the paper by Bahr K. 1991, Phys. Earth Planet Int., 66, 24-38. In this paper a systematic numerical experiment is given by using the corrected formulas by Bahr, when the disturbing body is considered as a near-surface electric dipole. Its position and orientation over a two-dimensional model is varied. The tensor decomposition carried out in this way perfectly recovers the phase values in the original 2D impedance tensor. 4.10p ON THE STABILITY OF MAGNETOTELLURIC TRANSFER FUNCTION ESTIMATES AND THE RELIABILITY OF THEIR VARIANCES Markus Eisel(1) and Gary D. Egbert(2) (1) GeoForschungsZentrum Potsdam, Now at S.E.S.A. AG, Germany (2) College of Oceanic and Atmospheric Sciences, Ocean Admin Bldg 104, Oregon State University, Corvallis, OR, 97331-5503, USA egbert@oce.orst.edu Using data from a continuously operating two-station magnetotelluric (MT) array in central California we study the stability of MT transfer function (TF) estimates and the reliability of error estimates obtained by several methods. Typical deviations of estimates based on a single day of data differ from the long term average TF by 2--3% for T<300s, increasing to about 10% for T=2000s. Day-to-day deviations from the average are largely random, though there is some evidence for small frequency independent variations in impedance amplitudes, suggestive of temporal variations in surface distortion. Comparison of estimated error bars to TF variability showed that where coherent noise effected TF estimates (T=10-100 s) the standard error bars for the robust estimator were too small by as much as a factor of two.Otherwise these error bars were consistent with the actual precision of the TF estimates. We also considered two approximate jackknife schemes. First, weights for the robust TF estimates were computed once with all data, and the jackknife was applied to the final weighted LS estimate. These error bars were larger, but still too small in the 10--100 second bias band and now too large at other periods. We also tried a subset deletion jackknife, applying the full robust procedure with subsets of data deleted. Provided large subsets (5%) were deleted this yielded significantly more realistic error bars in the bias band, but error bars at other periods were now often much too large. The jackknifed error bars were thus more conservative, though not necessarily more reliable. 4.11p GENERALIZED RICCATI EQUATION FOR 1-D MT IMPEDANCES IN ANISOTROPIC STRUCTURES S. Kovacikova and J. Pek Geophysical Institute, Academy of Sciences of the Czech Republic, Bocni II/1401, 14131 Praha 4, Czech Republic SVK@ig.cas.cz In the 1-D magnetotelluric theory, a Riccati equation governs the depth variations of the impedance, or admittance, for a given distribution of the electrical conductivity. This equation can be used to compute the surface MT functions for generally piecewise continuous conductivity profiles.In case of a simple layered medium, it provides the cassical recurrent formulae for recalculating the impedances between the individual layer boundaries. We present an extended version of the Riccati differential equations for generally anisotropic 1-D structures, both for the case of a plane wave incident field and a non-homogenous primary field. In the anisotropic case, all elements of the impedance tensor are present and, consequently, a system of four coupled Riccati equations has to be considered. By applying standard methods for the numerical solution of systems of ordinary differential equations to the Riccati system, we demonstrate that, besides providing the possibility to model 1-D anisotropic transition zones, the Riccati formulation of the direct problem has a better numerical perfofmance as compared to standard matrix methods for field components. A synthetic study on the influence of a depth- variable regional strike on MT decomposition results is presented, with the variable strike simulated by a variable anisotropy within the 1-D section. 4.12p A COMPLEX DISTORTION MATRIX TENSOR DECOMPOSITION, A NEW METHOD TO DETECT AND REMOVE LOCAL STRONG CURRENT CHANNELLING Pamela Lezaeta Freie Universitat Berlin, Germany pamela@geophysik.fu-berlin.de Telluric and magnetic galvanic effects due to local heterogeneities are expressed in a matrix of complex numbers. Each complex number can be expressed as function of the telluric and magnetic distortion parameters already defined in known decomposition schemes. The solution is direct, with no need of least square fitting to find the regional response and strike angle. Instead, the minimum standard deviation of each averaged distortion parameter is necessary condition to determined a frequency independent strike angle. The assumption of local strong current channelling affecting the regional fields, due to a high conductivity shallow elongated anomaly or vertical thin dike oriented in an arbitrary local azimuth, is also possible to be inferred by quantifying the "phase deviations" of the tensor elements. A frequency dependent regional strike angle can be obtained if this specific distortion assumption is valid. If both frequency independent and dependent strike angles are similar along the period band, then the assumption of strong current channelling under galvanic limits has been verified. The measured data can be corrected, i.e., the regional impedance tensor as well as the magnetic transfer function can be estimated therefrom. 4.13p AN APPROACH TO CORRECT GALVANIC DISTORTIONS EFFECTS ON 3D ENVIRONMENTS Juanjo Ledo(1), Pilar Queralt(2), and Anna Marti(2) (1) Geological Survey of Canada, 615 Booth St. K1A 0E9, Ottawa, Ontario, Canada (2) Dept. de Geodinamica i Geofisica, Universitat de Barcelona, Barcelona, Spain ledo@cg.nrcan.gc.ca Galvanic distortion produces a frequency-dependent effect on the measured data over a 3D structure with the result that even closed sites, responding to the same regional electrical structure, could shown significant different impedance tensor components. The measured impedance tensor over a 3D regional structure can be expressed as the result of the product between the distortion matrix and the regional 3D-impedance tensor. With some algebra the measured impedance components, 8 values, can be expressed as a function of the regional components ,8 unknowns, and the twist and shear components of the distortion matrix ,2 unknowns. To retrieve the regional 3D-impedance tensor from the measured one an approximation is needed in order to reduce the number of unknowns to eight. Another important aspect when dealing with complex data is to define the target we want to study given both the impossibility to model all the structures and to not over interpret the data. In this case we want to retrieve the regional-deeper structures. Taking into account the previous comments, the approach presented in this paper consist of supposing that the earth response can be modeled as a 2D structure for the shorter periods. Then the Groom-Bailey decomposition technique is applied to determine the twist and shear parameters, reducing the number the unknowns to eight. We present different examples, with synthetic and real data to study the effects of galvanic distortions over 3D and quasi-3D ,finite strike, regional structures and how to retrieve the regional response. 4.14p STATIC SHIF LEVELLING USING THE TIPPER Juanjo Ledo(1), Anna Gabas(2), and Alex Marcuello(2) (1) Geological Survey of Canada, 615 Booth St. K1A 0E9, Ottawa, Ontario, Canada (2) Dept. de Geodinamica i Geofisica, Universitat de Barcelona, Barcelona, Spain ledo@cg.nrcan.gc.ca The build up of current and charge distributions in the near surface inhomogeneities produces inductive and galvanic distortions respectively. The local current distribution affects both the regional electric and magnetic fields. However, the inductive effect rapidly vanishes with increasing period and these will not be taken into account in this work. The galvanic effect is due to electric field redistribution and there is no temporal phase variation between the regional electric field and the secondary electric field produced by the local, near-surface inhomogeneities. We assume in this work that the electric field generated by these charge distributions is curl free. Using the physical properties of the distorted electric field, we propose a method to estimate the value of the static shift from the measured vertical magnetic field through the tipper. Although this method is valid for 2D and 3D structures, we present here the 2D case. Two examples with synthetic data and one using real data show the utility of this method. 4.15p FREELY DECAYING CURRENT MODES: AN EXAMPLE FOR A COMPLETE SET OF EINGENFUNCTIONS AND POSSIBLE APPLICATIONS Peter Weidelt and Christiane Stuntebeck Institute of Geophysics and Meteorology, Technical University of Braunschweig, Braunschweig, Germany p.weidelt@tu-bs.de When the sources of a given current distribution in a conductor are switched off, the current freely decays and its variation in space and time can be represented as a superposition of ortho-normal free-decay modes. Each mode decays exponentially (with the decay constant as eigenvalue) and leaves the spatial structure unchanged with time. The striking feature of free-decay modes is their independence of any source: once the set of eigenfunctions is known for a given conductivity distribution, the response of the conductor to any excitation both in the frequency and time domain can be represented by a 'simple' superposition of the modes. This feature might make them attractive for computer-intensive modelling, required for instance in the interpretation of airborne data. In this contribution we are giving as an example the complete set of the vectorial free-decay modes for the simple model of a perfectly conducting half-plane in a uniform full-space. These modes are essentially the responses for incident plane waves. For specific sources we address the practically important convergence problem of the superposition if only a limited number of modes is available. 4.16p NUMERICAL MODELLING OF TRANSIENT-ELECTROMAGNETIC FIELDS IN THREE-DIMENSIONAL CONDUCTORS: A COMPARISON OF METHODS Peter Weidelt Institute of Geophysics and Meteorology, Technical University Braunschweig, Mendelssohnstr. 2-3, D- 38106 Braunschweig, Germany p.weidelt@tu-bs.de Advanced interpretative methods in applied geo-electromagnetism require the response of pulsed controlled sources in a three- dimensional conductivity structure. Despite recent advances in numerical modelling, this problem still presents a challenging task. In this contribution I assess the performance of two different algorithms, namely a stabilized time-stepping algorithm with artificial introduction of quasi-displacement currents (Dufort-Frankel method) and a spectral finite difference method based on the Lanczos algorithm (first proposed by Druskin & Knizhnerman, 1988). Both algorithms use staggered grids and confine the domain of computation to the conducting earth half-space by applying an integral boundary condition. Codes for these algorithms are not public domain and therefore have to be rewritten. These new versions are distinguished by a more efficient and accurate implementation of the surface boundary condition. The advantage of the time-stepping algorithm is its formal simplicity, which unfortunately has to be seen in relative terms because of the small time step required to conserve the diffusive nature of the solution in the presesence of displacement currents. Therefore late times require lengthy computing times and are difficult to achieve. The late times also pose a problem for the spectral difference approach, which require a great number of approximate free decay modes. The two programs will be freely distributed and thus will fill be of service for the induction community. 4.17p NEW ADVANCES IN A FAST 3D ELECTROMAGNETIC INTEGRAL EQUATION SOLUTION D.B. Avdeev(1), A.V. Kuvshinov(1), O.V. Pankratov(1), and G. Newman(2) (1) Geoelectromagnetic Research Institute, Russian Academy of Sciences, 142092 Troitsk, Moscow Region, Russia (2) Sandia National Laboratories, Org. 6116, P.O. Box 5800, Albuquerque, NM 87185-0750, USA d.avdeev@mtu-net.ru Electromagnetic (EM) modelling based on the modified iterative dissipative method (MIDM) has proven to be an effective way of solving the 3D scattering problems in geophysics. The method was originally presented by Singer (1995) for a quasi-static field in isotropic formations, and later extended to formations with displacement currents and anisotropy in (Pankratov et al., 1995, 1997; Singer & Fainberg, 1995, 1997). In this method, conventional scattering equation is modified to an integral equation (IE) which is then solved by simple iteration that is equivalent to the Neumann series (NS) summation. Since this modified IE has the property that the norm of the discrete integral operator is less than one, simple iteration monotonically converges for any frequency and for any contrast of electrical resistivity. Based on this method, a number of 3D numerical solutions have been recently developed for variety of EM applications (cf. Avdeev et al., 1997, 1998, 1999; Zhdanov & Fang, 1997; Koyama & Utada, 1998; Kuvshinov et al., 1999; Singer, Mezzatesta & Wang, 1999; Zhdanov et al., 1999). Using modern methods for solving linear systems, we have succeeded in further accelerating these fast IE solutions. Model studies demonstrate that an order of magnitude speed up is possible compared to conventional MIDM solutions. It is demonstrated also that with higher contrasting models the greater the acceleration that can be achieved. Our findings demonstrate this fact for contrasts as large as fifty thousand to one.