**Geophysics in Near-Surface Investigations**

## Jadwiga A. Jarzyna et al.\*

*AGH University of Science and Technology, Faculty of Geology Geophysics and Environmental Protection, Krakow, Poland* 

## **1. Introduction**

44 New Achievements in Geoscience

Salem, A., Ravat, D., Mushayandebvu, M. F., & Ushijima, K. (2004). Linearized least-squares

Silva, J. B. C. (1989). Transformation of nonlinear problems into linear ones applied to the

Silva, J. B. C., Oliveira, A. S., & Barbosa, V. C. F. (2010). Gravity inversion of 2D basement

Stanley, J. M. (1977) Simplified magnetic interpretation of the geologic contact and thin dike.

Tarantola, A. (2005). *Inverse problem theory and methods for model parameter estimation,* SIAM,

Thompson, D. T. (1982). EULDPH: A new technique for making computer-assisted depth

Tikhonov, A. N., & Arsenin, V. Y. (1977). *Solutions of ill-posed problems,* Winston & Sons,

Zhang, J., Zhong, B., Zhou, X., & Dai, Y. (2001). Gravity anomalies of 2-D bodies with variable density contrast. *Geophysics,* Vol. 66, No. 3, pp. 809–813, ISSN 0016-8033. Zhang, C., Mushayandebvu, M. F., Reid, A. B., Fairhead, J. D., & Odegard, M. E. (2000).

Zhdanov, M. S. (2002). *Geophysical inverse theory and regularization Problems,* Elsevier, ISBN 0-

*Geophysics,* Vol. 69, No. , pp. 783–788, ISSN 0016-8033.

*Geophysics,* Vol. 42, No. 6, pp. 1236-1240, ISSN 0016-8033.

ISBN 978-0-898715-72-9, Philadelphia, USA.

ISBN 0-521-43064-X, Washington, DC, USA.

512-520, ISSN 0016-8033.

444-51089-3, Amsterdam.

ISSN 0016-8033.

0016-8033.

method for interpretation of potential field data from sources of simple geometry.

magnetic field of a two-dimensional prism. *Geophysics,* Vol. 54, No. 1, pp. 114-121,

relief using entropic regularization. *Geophysics,* Vol. 75, No. 3, pp. I29–I35, ISSN

estimates from magnetic data. *Geophysics,* Vol. 47, No. 1, pp. 31-37, ISSN 0016-8033.

Euler deconvolution of gravity tensor gradient data. *Geophysics,* Vol. 65, No. 2, pp.

Environmental and engineering geophysics are new branches of geophysics that focus on data collection and an analysis of the present day environment. Noninvasive, nondestructive, and inexpensive geophysical technologies are now capable of providing detailed 3D and 4D representations of the subsurface. These advances in geophysical survey techniques have not only enhanced traditional applications, such as prospecting geophysics, but also made way for development of new methodologies in near surface investigations. Environmental and engineering geophysics includes mapping of shallow geological formations, prospecting for resources located in the near surface, monitoring aquifers, mapping anthropogenic effects on the environment, and surveying for civil engineering and archaeological projects. Environmental and engineering geophysics measures parameters such as salinity of subsurface fluids, levels of radioactivity, and other soil properties affected by industrial activity.

The increasing demand for sustainable development as well as community efforts to preserve and conserve the environment have fostered the development of engineering technologies that are more environmentally friendly than traditional geophysical tools. Newer techniques are often used to monitor and protect the environment, as well as evaluate geotechnical risks of various natural and human - made structures. The increasing number of these modern geophysical applications, and of professionals specializing in them, bear witness to their emerging significance.

## **2. Petrophysical background for geophysical investigations (Jadwiga A. Jarzyna)**

Petrophysics is a subdiscipline of geophysics that addresses the physical properties and other parameters of rock and other subsurface material. Applied geophysics combines the theoretical foundation of a given geophysical method with empirical understanding of the materials involved. In the case of petrophysics, theoretical parameters of rock and other subsurface materials (i.e. their physical properties) are used to interpret empirical

<sup>\*</sup> Jerzy Dec, Jerzy Karczewski, Sławomir Porzucek, Sylwia Tomecka-Suchoń,

Anna Wojas and Jerzy Ziętek

*AGH University of Science and Technology, Faculty of Geology Geophysics and Environmental Protection, Krakow, Poland* 

Geophysics in Near-Surface Investigations 47

properties, and then describes their use in environmental and engineering geophysical applications. The most basic petrophysical property of rock is that of density. **Bulk density** is defined as the ratio of a rock's mass to its volume. It is lower or equal to matrix density (mineralogical density, specific density), which is defined as the ratio of the mass of a rock in powdered form, to its volume. Bulk density of sedimentary rocks, δb, depends on the density of minerals (matrix material), δm; porosity, Ф; and density of the pore space media,

(1 ) ( ) *b m f m mf*

 

Fig. 2.1 Permeability vs. porosity of the Rotliegend sandstones; a) datasets from different

wells; b) datasets ordered according FZI (Jarzyna et al., 2009).

Bulk density is intuitively inversely proportional to porosity, but it is specifically a function of matrix and pore media density. For example, gaseous hydrocarbons in the pore space of a reservoir rock will significantly decrease its bulk density. Most rock-forming minerals have bulk densities ranging from 2200 to 3500 kg/m3. Ore bodies and associated minerals have average densities of 4000 to 8000 kg/m3. Sedimentary rocks are composed of matrix materials that range in average density from 2500 to 2900 kg/m3, of fluids ranging from 800 to 1240 kg/m3, and of gases that have average densities of less than 1000 kg/m3 (Kobranova, 1989). Given variation in the densities of sedimentary components, the range in bulk density can be quite broad from lignite (1000 - 1300 kg/m3); clay (1300 - 2300 kg/m3); sand and gravel (1400 - 2300 kg/m3); loam (1500 - 2200 kg/m3); sandstone (2000 - 2800 kg/m3); shale (2300 - 2800 kg/m3); limestone (2300 - 2900 kg/m3); and dolomite (2400 - 2900 kg/m3); (Schön, 2004). **Porosity** is an important property of sedimentary rocks that often enters into petrophysical considerations. Total porosity is defined as the volume of free space per given volume of rock. Effective porosity refers only to the volume of connected pores, fractures, fissures and vugs. Porosity can be modelled as spherical volumes of various radii (Kobranova, 1989) but is more reliably measured in the laboratory using rock samples or well log data. Porosity is related to texture, mineral composition, sedimentary fabric and a) b

  (2.1)

δf (2.1).

observations (i.e. quantities, subsurface extent etc.). For example, geoeletrical investigations use Archie's law (theoretical) to interpret voltage and current measurements (empirical) to gain insight into the effective porosity and resistivity of the materials under study.

Petrophysics provides empirical information on rocks from laboratory measurements and field investigations. Laboratory measurements offer a direct way to determine rock parameters, but the results are not always directly applicable to the real world. The relatively small size of samples studied in typical laboratory experiments, and the complexity of the natural environment bedevil integration of laboratory results not just in applied geophysics, but in nearly all Earth science disciplines. High costs of drilling and limited number of cores can also constrain the practicality of investigations that rely heavily on parameters gathered through laboratory experiments.

In spite of these constraints, petrophysics aims to consistently integrate information from laboratory and field measurements. In subsurface investigations for example, the accuracy and resolution of geological profiles depend heavily on parameters provided by field measurements. Field measurements such as resistivity, measured by surface geoelectric methods and resistivity well logs, in turn reflect both the heterogeneity of subsurface rock formations and operational characteristics of the device used to perform the measurement. Given the variability among these factors, statistical methods can be used to synthesize and interpret data obtained using different methodologies. Empirical formulas generated through petrophysical investigations can also be used to standardize, calibrate and scale geophysical devices.

In practice, petrophysics often provides geophysical models with estimates of subsurface parameters. Estimated parameters can include bulk density, resistivity, dielectric permittivity, velocity of elastic waves, and magnetic susceptibility. In reservoir characterization for example, petrophysics parameterizes porosity, permeability and filtration factors. Accurate and consistent estimation of these parameters requires understanding of mineral composition, porosity, and properties of the interstitial media, (i.e. salinity, temperature, etc.). In larger scale studies of anisotropy, petrophysics provides estimates of structure, texture and facies type for materials under investigation.

Petrophysical analysis begins with the assumptions that the subsurface is heterogeneous and may undergo change. Subsurface materials may vary in terms of their geometry, size, composition, structure and other properties. Petrophysical and petrochemical variables address heterogeneity in these properties by averaging them over unit volumes. Subsurface materials are also affected by temperature, pressure and the resonant frequencies of physical fields used to evaluate them. Near surface investigations are less dependent on subsurface pressure and temperature conditions than traditional (deeper) geophysical applications, but these factors still come into play. Lack of compaction or varying degrees of consolidation in near surface materials can result in differing resistivity and bulk density for rocks of identical lithology, but occurring at different depths. Physical fields are also affected by consolidation and compaction. Loose subsurface sands and gravels, as well as sandstones and mudstones in the aeration zone have anomalously high resistivity for example, due to the air content of their pore space.

Because most of the near surface is composed of sediments and sedimentary rock, this chapter focuses on these materials. The chapter first outlines the basics of petrophysical

observations (i.e. quantities, subsurface extent etc.). For example, geoeletrical investigations use Archie's law (theoretical) to interpret voltage and current measurements (empirical) to

Petrophysics provides empirical information on rocks from laboratory measurements and field investigations. Laboratory measurements offer a direct way to determine rock parameters, but the results are not always directly applicable to the real world. The relatively small size of samples studied in typical laboratory experiments, and the complexity of the natural environment bedevil integration of laboratory results not just in applied geophysics, but in nearly all Earth science disciplines. High costs of drilling and limited number of cores can also constrain the practicality of investigations that rely heavily

In spite of these constraints, petrophysics aims to consistently integrate information from laboratory and field measurements. In subsurface investigations for example, the accuracy and resolution of geological profiles depend heavily on parameters provided by field measurements. Field measurements such as resistivity, measured by surface geoelectric methods and resistivity well logs, in turn reflect both the heterogeneity of subsurface rock formations and operational characteristics of the device used to perform the measurement. Given the variability among these factors, statistical methods can be used to synthesize and interpret data obtained using different methodologies. Empirical formulas generated through petrophysical investigations can also be used to standardize, calibrate and scale

In practice, petrophysics often provides geophysical models with estimates of subsurface parameters. Estimated parameters can include bulk density, resistivity, dielectric permittivity, velocity of elastic waves, and magnetic susceptibility. In reservoir characterization for example, petrophysics parameterizes porosity, permeability and filtration factors. Accurate and consistent estimation of these parameters requires understanding of mineral composition, porosity, and properties of the interstitial media, (i.e. salinity, temperature, etc.). In larger scale studies of anisotropy, petrophysics provides

Petrophysical analysis begins with the assumptions that the subsurface is heterogeneous and may undergo change. Subsurface materials may vary in terms of their geometry, size, composition, structure and other properties. Petrophysical and petrochemical variables address heterogeneity in these properties by averaging them over unit volumes. Subsurface materials are also affected by temperature, pressure and the resonant frequencies of physical fields used to evaluate them. Near surface investigations are less dependent on subsurface pressure and temperature conditions than traditional (deeper) geophysical applications, but these factors still come into play. Lack of compaction or varying degrees of consolidation in near surface materials can result in differing resistivity and bulk density for rocks of identical lithology, but occurring at different depths. Physical fields are also affected by consolidation and compaction. Loose subsurface sands and gravels, as well as sandstones and mudstones in the aeration zone have anomalously high resistivity for example, due to

Because most of the near surface is composed of sediments and sedimentary rock, this chapter focuses on these materials. The chapter first outlines the basics of petrophysical

estimates of structure, texture and facies type for materials under investigation.

gain insight into the effective porosity and resistivity of the materials under study.

on parameters gathered through laboratory experiments.

geophysical devices.

the air content of their pore space.

properties, and then describes their use in environmental and engineering geophysical applications. The most basic petrophysical property of rock is that of density. **Bulk density** is defined as the ratio of a rock's mass to its volume. It is lower or equal to matrix density (mineralogical density, specific density), which is defined as the ratio of the mass of a rock in powdered form, to its volume. Bulk density of sedimentary rocks, δb, depends on the density of minerals (matrix material), δm; porosity, Ф; and density of the pore space media, δf (2.1).

$$
\delta \mathcal{S}\_b = (1 - \Phi)\mathcal{S}\_m + \Phi \,\delta\_f = \mathcal{S}\_m - (\mathcal{S}\_m - \mathcal{S}\_f)\Phi \tag{2.1}
$$

Bulk density is intuitively inversely proportional to porosity, but it is specifically a function of matrix and pore media density. For example, gaseous hydrocarbons in the pore space of a reservoir rock will significantly decrease its bulk density. Most rock-forming minerals have bulk densities ranging from 2200 to 3500 kg/m3. Ore bodies and associated minerals have average densities of 4000 to 8000 kg/m3. Sedimentary rocks are composed of matrix materials that range in average density from 2500 to 2900 kg/m3, of fluids ranging from 800 to 1240 kg/m3, and of gases that have average densities of less than 1000 kg/m3 (Kobranova, 1989). Given variation in the densities of sedimentary components, the range in bulk density can be quite broad from lignite (1000 - 1300 kg/m3); clay (1300 - 2300 kg/m3); sand and gravel (1400 - 2300 kg/m3); loam (1500 - 2200 kg/m3); sandstone (2000 - 2800 kg/m3); shale (2300 - 2800 kg/m3); limestone (2300 - 2900 kg/m3); and dolomite (2400 - 2900 kg/m3); (Schön, 2004). **Porosity** is an important property of sedimentary rocks that often enters into petrophysical considerations. Total porosity is defined as the volume of free space per given volume of rock. Effective porosity refers only to the volume of connected pores, fractures, fissures and vugs. Porosity can be modelled as spherical volumes of various radii (Kobranova, 1989) but is more reliably measured in the laboratory using rock samples or well log data. Porosity is related to texture, mineral composition, sedimentary fabric and

a) b

Fig. 2.1 Permeability vs. porosity of the Rotliegend sandstones; a) datasets from different wells; b) datasets ordered according FZI (Jarzyna et al., 2009).

Geophysics in Near-Surface Investigations 49

The specific resistivity of air equals 1014 Ohmm, and is similar to that of hydrocarbons. Electrolytes are ionic conductors. Conductivity of porous, fractured rocks is caused by ionic movement in subsurface fluids. Under these conditions, resistivity depends primarily on mineral type, rock fabric, temperature of the subsurface fluid, and the volume of connected pores. The resistivity of pure water (containing only H+ and OH- ions) is very high ( > 105 Ohms; Hearst et al., 2000). Conductivity, σw, and resistivity of subsurface water (ρw),

1

In equations 2.3 and 2.4, n is the number of components in the subsurface water, is the degree of dissociation, c is the concentration, z is the valence, and v is the mobility. Factors of α and v depend on temperature. Models often assume that pore fluid is a simple sodiumchloride solution. Formula 2.4 shows that pore fluid resistivity is temperature dependent. The variables t1, t2 are temperatures given in Celsius degrees. The relative dielectric permittivity of pore water, εr=81 (at t=21oC) differs from that of oil and gas (εr < 3). Salt concentration, cmol (given in moles per liter), exerts a relatively small influence on the relative dielectric permittivity of pore fluids. The relative dialectric permittivity of pore

\_ \_ 13 1.065 0.03006 *r r pure water mol mol mol*

Dielectric permittivity, εr, the velocity of electromagnetic waves, VEM, and EM attenuation, AEM, are used in georadar investigations. The dielectric permittivity of porous rocks strongly depends on the degree of water saturation. The salinity of pore water does not strongly affect εr (2.5), but does influence conductivity and attenuation of electromagnetic waves. Selected values for the electromagnetic properties of different pore media as well as surface

leum ice ice

sand (dry)

εr 25.8 56.2 2.0-2.7 2.0-2.2 3.1 4,15 3.2 1.2 1.55

εr 80 80 3-4 5-30 4 25 1 5-40 60-80 VEM [m/ns] 0.03 0.03 0.16 0.07 0.15 0.06 0.30 0.05-0.13 0.03-0.04 AEM [dB/m] 0.1 1000 0.01 1-100 0.01 0.03-0.3 0 1-300 0.3 Table 2.2 Relative dielectric permittivity of selected materials; \*from pure distilled water (-12ºC, 106 Hz), \*\*from pure distilled water (-12º C, 109 Hz), ^freshly fallen, hard packed, followed by rain (-20º C, 106 - 109 Hz), ^^freshly fallen, hard packed, followed by rain (-6° C, 106 Hz) (after Schön, 2004; # from http://www.physics.utoronto.ca/~exploration/courses).

silts #

 *czv*

1

2 21.5 () () 21.5 *w w <sup>t</sup> t t*

*t*

(2.3)

2 3

*cc c* (2.5)

\*

sand

ice \*\* snow ^

(wet) air clay peat

snow ^^

(2.4)

1 *<sup>n</sup> w iii i w i*

2 1

 

fluids, εr, is related to cmol as follows (Schön, 2004):

and subsurface materials are presented below in Table 2.2.

rine oil petro-

ice #

glyce-

water (saline)

 

Material etha-

Material water

nol

(fresh)

depends on several factors (2.3; 2.4):

degree of lithification. Clay content reduces porosity. Permeability is the ability of a material to transmit pore media, which depends heavily on the material's porosity. A functional alternative to calculating a rock's hydraulic properties as separate porosity and permeability variables is the Flow Zone Index (FZI). FZI depends on permeability and porosity, specifically considers tortuosity of pore space and specific surfaces without parameterizing either of these quantities. The advantages of using FZI over porosity and permeability are presented in figure 2.1 (Tiab & Donaldson, 2004; Jarzyna et al., 2009).

**Resistivity** and **dielectric permittivity** refer to electric properties of rocks. Both properties depend on mineral composition, porosity and properties of the pore space media. Dielectric permittivity is inversely proportional to the electromagnetic (EM) field frequencies often used in geophysical measurements. Factors such as water saturation may also influence the dielectric permittivity of subsurface materials. Frequency dependence makes it a dispersive parameter. Because the dielectric permittivity of water is relatively high compared to permittivity values for minerals and hydrocarbons, this parameter can be used to distinguish hydrocarbon plumes in water saturated materials. In terms of resistivity, silicates and carbonates have high specific resistivities (>109 Ohmm; Kobranova, 1989) and are thus classified as insulators. Sulfides and some oxides meanwhile are classified as mineral conductors (specific resistivities between 10-6 and 10-1 Ohmm), as are native metals (10-8 to 10-5 Ohmm). Variation in the specific resistivities of various minerals can be due to the presence of impurities, subsurface structural effects, and anisotropy. Relative dielectric permittivity, εr, is defined as the dielectric permittivity of rock, ε, divided by the dielectric permittivity of vacuum, ε0. For most rock forming minerals, εr ranges between 4 and 10 (Table 2.1). Relative dielectric permittivity is correlated to bulk density (2.2). The generalized formula (2.2) applies to most minerals with the exception of hydrated forms of montmorillonite and members of the sulfide group (Schön, 2004).


$$
\varepsilon\_r = \left(1.93 \pm 0.17\right)^{\delta\_b} \tag{2.2}
$$

Table 2.1 Petrophysical properties of minerals and rocks (\*after Halliburton, 1991; \*\* Schön, 2004; ^Kobranova, 1989; #calculated on the basis of averaged data).

degree of lithification. Clay content reduces porosity. Permeability is the ability of a material to transmit pore media, which depends heavily on the material's porosity. A functional alternative to calculating a rock's hydraulic properties as separate porosity and permeability variables is the Flow Zone Index (FZI). FZI depends on permeability and porosity, specifically considers tortuosity of pore space and specific surfaces without parameterizing either of these quantities. The advantages of using FZI over porosity and permeability are

**Resistivity** and **dielectric permittivity** refer to electric properties of rocks. Both properties depend on mineral composition, porosity and properties of the pore space media. Dielectric permittivity is inversely proportional to the electromagnetic (EM) field frequencies often used in geophysical measurements. Factors such as water saturation may also influence the dielectric permittivity of subsurface materials. Frequency dependence makes it a dispersive parameter. Because the dielectric permittivity of water is relatively high compared to permittivity values for minerals and hydrocarbons, this parameter can be used to distinguish hydrocarbon plumes in water saturated materials. In terms of resistivity, silicates and carbonates have high specific resistivities (>109 Ohmm; Kobranova, 1989) and are thus classified as insulators. Sulfides and some oxides meanwhile are classified as mineral conductors (specific resistivities between 10-6 and 10-1 Ohmm), as are native metals (10-8 to 10-5 Ohmm). Variation in the specific resistivities of various minerals can be due to the presence of impurities, subsurface structural effects, and anisotropy. Relative dielectric permittivity, εr, is defined as the dielectric permittivity of rock, ε, divided by the dielectric permittivity of vacuum, ε0. For most rock forming minerals, εr ranges between 4 and 10 (Table 2.1). Relative dielectric permittivity is correlated to bulk density (2.2). The generalized formula (2.2) applies to most minerals with the exception of hydrated forms of

(1.93 0.17) *<sup>b</sup>*

Quartz 5.492 4.119 2650 Sandstone^ 0.8-4.5 30-100 0.06-0.6 Calcite 6.403 3.436 2710 6.35; 7.5-8.7 Sandstone# 1.4-4.3 3-41 0.13-0.33 Orthoclase 5.690 3.260 2570 5.6; 4.5-6.2 Sandy shale^ 1.45-5.18 5-69 0.12-0.21 Dolomite 7.007 4.293 2870 7.46; 6.3-8.2 Sand (dry)# 0.2-1.0 0.03-0.72 0.405 Anhydrite 6.096 3.126 2960 6.5; 5.7-6.7 Sand (wet)# 0.8-2.2 0.55-4.18 0.405 Siderite 6.959 3.590 3960 9.3; 5.2-7.4 Clay^ 0.3-3 30 0.25-0.45 Pyrite 8.021 5.166 5020 33.7-81 Clay# 1.0-2.5 0.78-4.91 0.405 Hematite 6.626 4.233 5270 25 Shale (slate)^ 2.3-6.65 24-72 0.17 Magnetite 4.175 1.966 5180 33.7-81 Limestone^ 1.0-7.0 13-175 0.18-0.31 Kaolinite 1.438 0.929 2610 11.8; 9.1 Limestone# 5.9-6.1 55-63 0.34-0.354 Biotite 6.220 3.717 3010 6.3; 6.2-9.3 Marl^ 1.3-4.5 10-135 0.11-0.23 Halite 4.549 2.628 2170 5.9; 5.7-6.2 Granite# 5.5-5.9 56-64 0.325 Table 2.1 Petrophysical properties of minerals and rocks (\*after Halliburton, 1991; \*\* Schön,

[MHz] Rock VP

(2.2)

[km/s]

E [GPa] <sup>σ</sup>

presented in figure 2.1 (Tiab & Donaldson, 2004; Jarzyna et al., 2009).

montmorillonite and members of the sulfide group (Schön, 2004).

\*δ<sup>b</sup> [kg/m3]

2004; ^Kobranova, 1989; #calculated on the basis of averaged data).

Mineral \*VP

[km/s]

\*VS [km/s] *r*

\*\*ε<sup>r</sup>

The specific resistivity of air equals 1014 Ohmm, and is similar to that of hydrocarbons. Electrolytes are ionic conductors. Conductivity of porous, fractured rocks is caused by ionic movement in subsurface fluids. Under these conditions, resistivity depends primarily on mineral type, rock fabric, temperature of the subsurface fluid, and the volume of connected pores. The resistivity of pure water (containing only H+ and OH ions) is very high ( > 105 Ohms; Hearst et al., 2000). Conductivity, σw, and resistivity of subsurface water (ρw), depends on several factors (2.3; 2.4):

$$
\sigma\_w = \frac{1}{\rho\_w} \equiv \sum\_{i=1}^n \alpha\_i c\_i z\_i v\_i \tag{2.3}
$$

$$
\rho\_w(t\_2) = \rho\_w(t\_1) \frac{t\_1 + 21.5}{t\_2 + 21.5} \tag{2.4}
$$

In equations 2.3 and 2.4, n is the number of components in the subsurface water, is the degree of dissociation, c is the concentration, z is the valence, and v is the mobility. Factors of α and v depend on temperature. Models often assume that pore fluid is a simple sodiumchloride solution. Formula 2.4 shows that pore fluid resistivity is temperature dependent. The variables t1, t2 are temperatures given in Celsius degrees. The relative dielectric permittivity of pore water, εr=81 (at t=21oC) differs from that of oil and gas (εr < 3). Salt concentration, cmol (given in moles per liter), exerts a relatively small influence on the relative dielectric permittivity of pore fluids. The relative dialectric permittivity of pore fluids, εr, is related to cmol as follows (Schön, 2004):

$$
\boldsymbol{\varepsilon}\_{r} = \boldsymbol{\varepsilon}\_{r\_{-}\boldsymbol{\mu}ure\\_water} - 13 \cdot \boldsymbol{\varepsilon}\_{mol} + 1.065 \cdot \boldsymbol{\varepsilon}\_{mol}^2 - 0.09006 \cdot \boldsymbol{\varepsilon}\_{mol}^3 \tag{2.5}
$$

Dielectric permittivity, εr, the velocity of electromagnetic waves, VEM, and EM attenuation, AEM, are used in georadar investigations. The dielectric permittivity of porous rocks strongly depends on the degree of water saturation. The salinity of pore water does not strongly affect εr (2.5), but does influence conductivity and attenuation of electromagnetic waves. Selected values for the electromagnetic properties of different pore media as well as surface and subsurface materials are presented below in Table 2.2.


Table 2.2 Relative dielectric permittivity of selected materials; \*from pure distilled water (-12ºC, 106 Hz), \*\*from pure distilled water (-12º C, 109 Hz), ^freshly fallen, hard packed, followed by rain (-20º C, 106 - 109 Hz), ^^freshly fallen, hard packed, followed by rain (-6° C, 106 Hz) (after Schön, 2004; # from http://www.physics.utoronto.ca/~exploration/courses).

Geophysics in Near-Surface Investigations 51

Fig. 2.3 Velocity and bulk density vs. porosity; Miocene sandstones and anhydrites;

**Velocities of elastic waves** used in seismic surveys (longitudinal P-waves, and shear S-waves) strongly depend on porosity and elastic moduli of the subsurface matrix, and on the type and properties of the pore space media. P- and S-wave velocities, along with bulk density and lithological information are used to calibrate seismic profiles. The velocity ratio, VP/VS, and acoustic impedance can indicate differences in lithology. Many petrophysical parameters are strongly depth dependent, and wave velocities in particular can be obscured by depth dependent changes in materials (Fig. 2.2). Porosity typically reduces velocity (Fig. 2.3.), as described by the Wyllie equation, the Raymer-Hunt-Gardner formula, and by others (Schön, 2004). A velocity correction for lack of compaction is necessary for empirical analysis of unconsolidated near-surface materials. Clays lower the elastic moduli of rocks and thus significantly reduce velocities of elastic waves. Velocities in clay-rich materials are lower than those in quartz- or calcite- rich materials (Table 2.1) while the compressibility of clay is higher than that of quartz and calcite. Clay properties are also strongly influenced by water volume. The above factors lead to complex relationships among the physical properties of subsurface materials that

Wave velocities in frozen unconsolidated sediments (permafrost) are similar to velocities in consolidated sediments, as ice can serve as a cement. Pronounced differences in the VP/VS ratio for gas- and water-saturated rocks are useful in energy exploration and cases of subsurface hydrocarbon contamination. Differences in the VP/VS ratio for gas and liquid saturated materials arise from the differential influence of pore space media on P- and S-wave velocities (differing effects on VP and VS, respectively). Velocities (for Poisson ratio,

O.F.- units older than Permian.

have high clay and water content.

For rocks that have no clay minerals, Archie's laws (2.6 and 2.7) describe the relationships among the resistivity of fully saturated rock (R0), the resistivity of rock partially saturated with water and hydrocarbons (Rt), effective porosity (Фeff), resistivity of formation water (Rw), and water saturation (Sw) or hydrocarbon saturation (1-Sw). Relationship I vs. Sw (2.7; Archie law) may be nonlinear on a bi-logarithmic scale, especially if there are several different brine or porosity systems within the subsurface formation, or if the formation has a significant fine-grained lithological component (≤ mud sized particles).

$$F = \frac{R\_0}{R\_w} = \frac{1}{\Phi\_{ef}{}^m} \tag{2.6}$$

$$I = \frac{R\_t}{R\_0} = \frac{1}{S\_w \, ^\circ} \tag{2.7}$$

In equations 2.6 and 2.7, F is the formation factor, I is the resistivity index, m is a cementation factor ranging from 1.3 for unconsolidated sediments, to 2.2 for lithified, nonporous, non-permeable rock, and n is the saturation exponent, ranging from 1.12 to 2.55 for sandstone, and from 1.1 to 2.38 for limestone. The variables m and n are empirically determined (Tiab & Donaldson, 2004).

Fig. 2.2 Bulk density vs. depth; sandstone, S, limestone, L, anhydrite, A, Miocene N-E north-eastern Carpathian Foredeep.

For rocks that have no clay minerals, Archie's laws (2.6 and 2.7) describe the relationships among the resistivity of fully saturated rock (R0), the resistivity of rock partially saturated with water and hydrocarbons (Rt), effective porosity (Фeff), resistivity of formation water (Rw), and water saturation (Sw) or hydrocarbon saturation (1-Sw). Relationship I vs. Sw (2.7; Archie law) may be nonlinear on a bi-logarithmic scale, especially if there are several different brine or porosity systems within the subsurface formation, or if the formation has a

> <sup>0</sup> 1 *<sup>m</sup> <sup>w</sup> ef*

0 1 *<sup>t</sup> n w*

In equations 2.6 and 2.7, F is the formation factor, I is the resistivity index, m is a cementation factor ranging from 1.3 for unconsolidated sediments, to 2.2 for lithified, nonporous, non-permeable rock, and n is the saturation exponent, ranging from 1.12 to 2.55 for sandstone, and from 1.1 to 2.38 for limestone. The variables m and n are empirically

*<sup>R</sup>* (2.6)

*<sup>R</sup> <sup>S</sup>* (2.7)

*<sup>R</sup> <sup>F</sup>*

*<sup>R</sup> <sup>I</sup>*

Fig. 2.2 Bulk density vs. depth; sandstone, S, limestone, L, anhydrite, A, Miocene N-E -

significant fine-grained lithological component (≤ mud sized particles).

determined (Tiab & Donaldson, 2004).

north-eastern Carpathian Foredeep.

Fig. 2.3 Velocity and bulk density vs. porosity; Miocene sandstones and anhydrites; O.F.- units older than Permian.

**Velocities of elastic waves** used in seismic surveys (longitudinal P-waves, and shear S-waves) strongly depend on porosity and elastic moduli of the subsurface matrix, and on the type and properties of the pore space media. P- and S-wave velocities, along with bulk density and lithological information are used to calibrate seismic profiles. The velocity ratio, VP/VS, and acoustic impedance can indicate differences in lithology. Many petrophysical parameters are strongly depth dependent, and wave velocities in particular can be obscured by depth dependent changes in materials (Fig. 2.2). Porosity typically reduces velocity (Fig. 2.3.), as described by the Wyllie equation, the Raymer-Hunt-Gardner formula, and by others (Schön, 2004). A velocity correction for lack of compaction is necessary for empirical analysis of unconsolidated near-surface materials. Clays lower the elastic moduli of rocks and thus significantly reduce velocities of elastic waves. Velocities in clay-rich materials are lower than those in quartz- or calcite- rich materials (Table 2.1) while the compressibility of clay is higher than that of quartz and calcite. Clay properties are also strongly influenced by water volume. The above factors lead to complex relationships among the physical properties of subsurface materials that have high clay and water content.

Wave velocities in frozen unconsolidated sediments (permafrost) are similar to velocities in consolidated sediments, as ice can serve as a cement. Pronounced differences in the VP/VS ratio for gas- and water-saturated rocks are useful in energy exploration and cases of subsurface hydrocarbon contamination. Differences in the VP/VS ratio for gas and liquid saturated materials arise from the differential influence of pore space media on P- and S-wave velocities (differing effects on VP and VS, respectively). Velocities (for Poisson ratio,

Geophysics in Near-Surface Investigations 53

Iron-oxide minerals are often used as paleoenvironmental proxies due to their high magnetic susceptibility, and because environmental conditions can be interpreted from iron

**3. Near surface georadar investigations (Jerzy Karczewski & Jerzy Ziętek)** 

In recent years, ground penetrating radar (GPR) has become the most popular geophysical method for near surface investigations. GPR belongs to a group of radio wave methods which evaluate electromagnetic wave propagation within a geological medium. The most popular GPR apparatus is the impulse type, which uses transmitting antennas to emit short electromagnetic pulses of 0.5 to 10 ns in the 10MHz to 6GHz frequency range (f). Under favorable conditions, low frequency antennas (f < 50 MHz) can be used to map the structure of the subsurface down to depths of a few dozen meters. Higher frequency pulses increase the resolution of shallow features, but decrease the depth range of the survey. EM waves from GPR devices can be attenuated, reflected or refracted. Reflection or refraction coefficients of EM waves for a given subsurface boundary layer will depend on different dielectric permittivities (ε) of the materials involved. Dielectric permittivity also influences velocities of the EM waves. Wave attenuation depends on wave frequency and the conductivity of the media, which in turn is strongly influenced by depth. While GPR devices are cheap and relatively easy to operate, GPR echograms are difficult to interpret and require knowledge of subsurface structure and petrophysical properties. Methods for processing GPR field data can enhance signal to noise ratios and the accuracy of reflector images. The most common GPR processing procedures include DC-Bias, time gain, filtration in the frequency domain, f/k filtration, deconvolution, migration and others (Annan, 2001;

Landslides can damage infrastructure and pose a variety of other environmental risks. Mass wasting during a landslide event strongly depends on topography, hydrology, the structure of underlying bedrock, soil and bedrock type, and other factors. Strategies for landslide risk mitigation require a detailed understanding of the internal structure of the landslides, especially the slide surface. GPR surveys performed to image a landslide surface should profile the main landslide axis, as well as an axis perpendicular to the slide if possible.

For landslides that cut through inclined layers, a topographic correction is necessary to properly locate and map the slide surface. Stacking of multiple signals during data acquisition and maximizing the number stacked signals improves the significance and clarity of the echograms. If the terrain permits, several parallel profiles can be recorded to render a 3D map of the main slide surface. A small landside that destroyed a relatively new house in the south of Poland provides a case study of the methods described above. GPR was used to locate and map the landslide surface, which developed following a period of heavy rain in 1997 (Fig. 3.1). No landslide activity had been observed in the area prior to 1997, leading investigators to suspect that construction of the house and it adjacent swimming pool in 1995 may have disturbed the geotechnical balance of the property. Later drilling revealed the lithological composition of the subsurface rocks, and confirmed the

oxide composition and grain size distribution.

Daniels, 2004; Jol, 2009).

presence of slide surface.

**3.1 Landslide investigations** 

σ) and bulk density (for Young modulus, E, and others) can be used to calculate the dynamic elastic moduli of rocks. Sonic well logs record full, acoustic waveforms for *in situ* elastic wave velocity analysis. In unconsolidated, near-surface sediments, shear wave velocity can be lower than the velocity of borehole mud. Under these conditions, it is necessary to use a formula derived from the Stoneley wave equation in order to calculate S-wave velocity (Jarzyna et al., 2010). Sedimentary rocks are often anisotropic and well layered. When measured parallel to sedimentary layers, velocity, VP//, and Young's modulus, E//, are greater than VP┴ and E┴, measured perpendicular to layers. The anisotropic coefficient of longitudinal wave velocity can reach values of 1.2 to 1.3. These values can exceed transverse wave coefficients and anisotropic coefficients pertaining to Young's modulus (~1.1 to 2.0; Kobranova, 1989). In finely laminated rocks, such as shale-rich sandstones, elastic waves having wavelengths longer than the thickness of the individual laminae must be normalized in order to calculate the effective properties of the material.

**Magnetic susceptibility** is an important parameter for soils and shallow subsurface rocks containing ferrimagnetic minerals. The magnetic susceptibility of soil and rock depends on shape, size and concentration of susceptible minerals. Most minerals are dia- and paramagnetic. Antiferrimagnetic minerals are relatively rare and ferromagnetic minerals are very rare (Kobranova, 1989). Specific volumetric susceptibility, κ, and mass susceptibility, χ, respectively, are proportionality coefficients in the presented below mutual relationships:

$$
\vec{J}\_v = \mathbf{\bar{x}} \cdot \vec{H} \tag{2.8}
$$

$$
\vec{J}\_m = \mathcal{X} \cdot \vec{H} \tag{2.9}
$$

$$\mathcal{X} = \mathbf{x} \nmid \mathcal{S}\_m \tag{2.10}$$

In equations 2.8 and 2.9, *H* is magnetic field, *vJ* is specific volume, and *Jm* is specific mass magnetization.

The magnetic susceptibility of dia- and paramagnetic native elements is governed by their chemical properties. For multi-element minerals, magnetic susceptibility is governed by stoichiometric composition and the arrangement of elements in the crystal lattice. The specific volumetric magnetic susceptibility of diamagnetic minerals is negligible and negative in value, ranging from -5.024 to -1.63 x 10-4 SI units (Kobranova, 1989). The most common minerals found in sedimentary rocks are diamagnetic (i.e. quartz, calcite, feldspar, dolomite, gypsum, halite, etc). Trace minerals such as biotite, pyrite, ilmenite, siderite, chlorite, and clays are paramagnetic or paraferrimagnetic. Accessory minerals such as magnetite, maghemite, hematite, and goethite are ferro- or ferrimagnetic. Magnetite is one of the most susceptible minerals, with 1.25 to 25 or more SI units. Water and oil are diamagnetic (κw=χw= -0.9 10-5 and κo=χo= -1.04 10-5 SI units, respectively). The influence of the mineral matrix on the susceptibility of the pore fluid is insignificant. Gases, including gaseous hydrocarbons have much lower magnetic susceptibility than liquids. Oxygen is paramagnetic with a relatively low, κoxyg = 0.17 x 10-5 SI units. The magnetic susceptibility of air, κair, equals 0.04 x 10-5 SI units. The ferromagnetic properties of magnetized rocks are the basis of magnetic surveys, paleogeographic reconstructions, and the paleomagnetic record.

σ) and bulk density (for Young modulus, E, and others) can be used to calculate the dynamic elastic moduli of rocks. Sonic well logs record full, acoustic waveforms for *in situ* elastic wave velocity analysis. In unconsolidated, near-surface sediments, shear wave velocity can be lower than the velocity of borehole mud. Under these conditions, it is necessary to use a formula derived from the Stoneley wave equation in order to calculate S-wave velocity (Jarzyna et al., 2010). Sedimentary rocks are often anisotropic and well layered. When measured parallel to sedimentary layers, velocity, VP//, and Young's modulus, E//, are greater than VP┴ and E┴, measured perpendicular to layers. The anisotropic coefficient of longitudinal wave velocity can reach values of 1.2 to 1.3. These values can exceed transverse wave coefficients and anisotropic coefficients pertaining to Young's modulus (~1.1 to 2.0; Kobranova, 1989). In finely laminated rocks, such as shale-rich sandstones, elastic waves having wavelengths longer than the thickness of the individual laminae must be normalized

**Magnetic susceptibility** is an important parameter for soils and shallow subsurface rocks containing ferrimagnetic minerals. The magnetic susceptibility of soil and rock depends on shape, size and concentration of susceptible minerals. Most minerals are dia- and paramagnetic. Antiferrimagnetic minerals are relatively rare and ferromagnetic minerals are very rare (Kobranova, 1989). Specific volumetric susceptibility, κ, and mass susceptibility, χ, respectively, are proportionality coefficients in the presented below

> *vJ H*

*J H <sup>m</sup>* 

The magnetic susceptibility of dia- and paramagnetic native elements is governed by their chemical properties. For multi-element minerals, magnetic susceptibility is governed by stoichiometric composition and the arrangement of elements in the crystal lattice. The specific volumetric magnetic susceptibility of diamagnetic minerals is negligible and negative in value, ranging from -5.024 to -1.63 x 10-4 SI units (Kobranova, 1989). The most common minerals found in sedimentary rocks are diamagnetic (i.e. quartz, calcite, feldspar, dolomite, gypsum, halite, etc). Trace minerals such as biotite, pyrite, ilmenite, siderite, chlorite, and clays are paramagnetic or paraferrimagnetic. Accessory minerals such as magnetite, maghemite, hematite, and goethite are ferro- or ferrimagnetic. Magnetite is one of the most susceptible minerals, with 1.25 to 25 or more SI units. Water and oil are diamagnetic (κw=χw= -0.9 10-5 and κo=χo= -1.04 10-5 SI units, respectively). The influence of the mineral matrix on the susceptibility of the pore fluid is insignificant. Gases, including gaseous hydrocarbons have much lower magnetic susceptibility than liquids. Oxygen is paramagnetic with a relatively low, κoxyg = 0.17 x 10-5 SI units. The magnetic susceptibility of air, κair, equals 0.04 x 10-5 SI units. The ferromagnetic properties of magnetized rocks are the basis of magnetic surveys, paleogeographic reconstructions, and the paleomagnetic record.

 

is magnetic field, *vJ*

(2.8)

(2.9)

is specific volume, and *Jm*

/ *<sup>m</sup>* (2.10)

is specific mass

in order to calculate the effective properties of the material.

mutual relationships:

In equations 2.8 and 2.9, *H*

magnetization.

Iron-oxide minerals are often used as paleoenvironmental proxies due to their high magnetic susceptibility, and because environmental conditions can be interpreted from iron oxide composition and grain size distribution.

## **3. Near surface georadar investigations (Jerzy Karczewski & Jerzy Ziętek)**

In recent years, ground penetrating radar (GPR) has become the most popular geophysical method for near surface investigations. GPR belongs to a group of radio wave methods which evaluate electromagnetic wave propagation within a geological medium. The most popular GPR apparatus is the impulse type, which uses transmitting antennas to emit short electromagnetic pulses of 0.5 to 10 ns in the 10MHz to 6GHz frequency range (f). Under favorable conditions, low frequency antennas (f < 50 MHz) can be used to map the structure of the subsurface down to depths of a few dozen meters. Higher frequency pulses increase the resolution of shallow features, but decrease the depth range of the survey. EM waves from GPR devices can be attenuated, reflected or refracted. Reflection or refraction coefficients of EM waves for a given subsurface boundary layer will depend on different dielectric permittivities (ε) of the materials involved. Dielectric permittivity also influences velocities of the EM waves. Wave attenuation depends on wave frequency and the conductivity of the media, which in turn is strongly influenced by depth. While GPR devices are cheap and relatively easy to operate, GPR echograms are difficult to interpret and require knowledge of subsurface structure and petrophysical properties. Methods for processing GPR field data can enhance signal to noise ratios and the accuracy of reflector images. The most common GPR processing procedures include DC-Bias, time gain, filtration in the frequency domain, f/k filtration, deconvolution, migration and others (Annan, 2001; Daniels, 2004; Jol, 2009).

#### **3.1 Landslide investigations**

Landslides can damage infrastructure and pose a variety of other environmental risks. Mass wasting during a landslide event strongly depends on topography, hydrology, the structure of underlying bedrock, soil and bedrock type, and other factors. Strategies for landslide risk mitigation require a detailed understanding of the internal structure of the landslides, especially the slide surface. GPR surveys performed to image a landslide surface should profile the main landslide axis, as well as an axis perpendicular to the slide if possible.

For landslides that cut through inclined layers, a topographic correction is necessary to properly locate and map the slide surface. Stacking of multiple signals during data acquisition and maximizing the number stacked signals improves the significance and clarity of the echograms. If the terrain permits, several parallel profiles can be recorded to render a 3D map of the main slide surface. A small landside that destroyed a relatively new house in the south of Poland provides a case study of the methods described above. GPR was used to locate and map the landslide surface, which developed following a period of heavy rain in 1997 (Fig. 3.1). No landslide activity had been observed in the area prior to 1997, leading investigators to suspect that construction of the house and it adjacent swimming pool in 1995 may have disturbed the geotechnical balance of the property. Later drilling revealed the lithological composition of the subsurface rocks, and confirmed the presence of slide surface.

Geophysics in Near-Surface Investigations 55

planned according to known subsurface patterns, or the layout of the overlying structure. GPR surveys of the churches of Saint Peter and Saint Paul in Krakow were planned using the surface structure as a template (Fig. 3.3). In this example, the survey used shielded 250 MHz antennas and followed an aisle in the church that lay above known and accessible crypts. The longitudinal axes of the last two crypts (Fig. 3.3) are parallel to the profile. The echogram from the aisle was then used to interpret depth profiles of the nave of the church

Fig. 3.2 GPR of Saint Margaret's collegiate in central Poland; a) echogram, b) ruins of the

River embankments, levies and dams that secure reservoirs, holding ponds and other bodies of water, require periodic maintenance and upgrade. Even robustly constructed embankments degrade over time, usually from infiltration at the embankment's base. Water filtration across the body of the embankment results in soil suffusion and liquefaction.

Animal burrows, vegetation, and other biological modifications can also degrade embankments. Assessment of embankments is usually performed using geotechnical methods. These include drilling wells, conducting well tests, testing of well plugs in the laboratory, and *in situ* ground testing. Geomorphological factors affecting embankments can be ascertained via larger scale surveys of the area, using maps, photointerpretation,

Romanesque basement. RAMAC/GPR apparatus, 800 MHz antennas.

**3.3 GPR for controlling embankments** 

photogrammetry.

where unknown and inaccessible crypts were suspected to exist.

Fig. 3.1 GPR echogram of a landslide. RAMAC/GPR apparatus, 200 MHz antennas.

## **3.2 GPR in archeology**

More and more archeological investigations utilize geophysical methods of magnetometry, GPR, and geoelectrics. Features such as old moats, ramparts, walls, basements, cemeteries, residual strongholds, foundations, and tumuli have all been located using GPR methods. Archeological objects can be successfully identified using medium frequency antennas, as they are buried at depths of less than 10 m. GPR is non-invasive and thus offers the advantages of speed and precision in archeology, as well as the possibility of object recognition without digging an exploratory trench. Archeological interpretations can be optimized by doing GPR profiles of an entire grid system if terrain conditions permit. Under these conditions a GPR survey can significantly reduce the time and effort necessary to locate subsurface objects, thus lowering the costs of excavation. GPR can also help optimize design of research trenches and excavation plans for specific objects. Higher frequency GPR records can help differentiate small artifacts from other, less significant fragments of rock in the subsurface.

GPR is an especially effective, low impact method for surveys carried out inside historic buildings. Void space such as crypts, basements and tunnels can be detected as a large increase in EM wave velocity. As a result, the boundaries of void spaces are visible in echograms at apparent depths that are below the actual boundaries or floors of these features. Voids appear in echograms as low frequency impulses zones. Shielded antennas are often used for GPR measurements inside of buildings to reduce feedback from walls, ceilings and other objects (i.e. reflex anomalies), but unshielded antennas can also provide adequate results. Echograms collected inside buildings can contain a significant amount of noise from reflex anomalies, but these datasets have nevertheless facilitated numerous archeological discoveries.

GPR surveys of Saint Margaret's collegiate in central Poland for example, were used to determine whether a Romanesque basement excavated outside of the church continued beneath the floor of the historic structure above. A GPR survey using shielded antennas (800 MHz) confirmed that the basement extended beneath the church (Fig. 3.2). Due to the architectural continuity of many historical structures, archeological GPR profiles can be

Fig. 3.1 GPR echogram of a landslide. RAMAC/GPR apparatus, 200 MHz antennas.

More and more archeological investigations utilize geophysical methods of magnetometry, GPR, and geoelectrics. Features such as old moats, ramparts, walls, basements, cemeteries, residual strongholds, foundations, and tumuli have all been located using GPR methods. Archeological objects can be successfully identified using medium frequency antennas, as they are buried at depths of less than 10 m. GPR is non-invasive and thus offers the advantages of speed and precision in archeology, as well as the possibility of object recognition without digging an exploratory trench. Archeological interpretations can be optimized by doing GPR profiles of an entire grid system if terrain conditions permit. Under these conditions a GPR survey can significantly reduce the time and effort necessary to locate subsurface objects, thus lowering the costs of excavation. GPR can also help optimize design of research trenches and excavation plans for specific objects. Higher frequency GPR records can help differentiate small artifacts from other, less significant fragments of rock in

GPR is an especially effective, low impact method for surveys carried out inside historic buildings. Void space such as crypts, basements and tunnels can be detected as a large increase in EM wave velocity. As a result, the boundaries of void spaces are visible in echograms at apparent depths that are below the actual boundaries or floors of these features. Voids appear in echograms as low frequency impulses zones. Shielded antennas are often used for GPR measurements inside of buildings to reduce feedback from walls, ceilings and other objects (i.e. reflex anomalies), but unshielded antennas can also provide adequate results. Echograms collected inside buildings can contain a significant amount of noise from reflex anomalies, but these datasets have nevertheless facilitated numerous

GPR surveys of Saint Margaret's collegiate in central Poland for example, were used to determine whether a Romanesque basement excavated outside of the church continued beneath the floor of the historic structure above. A GPR survey using shielded antennas (800 MHz) confirmed that the basement extended beneath the church (Fig. 3.2). Due to the architectural continuity of many historical structures, archeological GPR profiles can be

**3.2 GPR in archeology** 

the subsurface.

archeological discoveries.

planned according to known subsurface patterns, or the layout of the overlying structure. GPR surveys of the churches of Saint Peter and Saint Paul in Krakow were planned using the surface structure as a template (Fig. 3.3). In this example, the survey used shielded 250 MHz antennas and followed an aisle in the church that lay above known and accessible crypts. The longitudinal axes of the last two crypts (Fig. 3.3) are parallel to the profile. The echogram from the aisle was then used to interpret depth profiles of the nave of the church where unknown and inaccessible crypts were suspected to exist.

Fig. 3.2 GPR of Saint Margaret's collegiate in central Poland; a) echogram, b) ruins of the Romanesque basement. RAMAC/GPR apparatus, 800 MHz antennas.

## **3.3 GPR for controlling embankments**

River embankments, levies and dams that secure reservoirs, holding ponds and other bodies of water, require periodic maintenance and upgrade. Even robustly constructed embankments degrade over time, usually from infiltration at the embankment's base. Water filtration across the body of the embankment results in soil suffusion and liquefaction.

Animal burrows, vegetation, and other biological modifications can also degrade embankments. Assessment of embankments is usually performed using geotechnical methods. These include drilling wells, conducting well tests, testing of well plugs in the laboratory, and *in situ* ground testing. Geomorphological factors affecting embankments can be ascertained via larger scale surveys of the area, using maps, photointerpretation, photogrammetry.

Geophysics in Near-Surface Investigations 57

way, GPR can optimize the more expensive assessment techniques, limiting their use to zones where GPR anomalies are apparent. The echogram shown in Fig. 3.4 was recorded on the Odra river embankment using 200 MHz unshielded antennas. The height of the embankment along the section analyzed was about 4 m. The goal of the investigation was to identify heterogeneities within the embankment and its basal materials. No meaningful heterogeneities were observed within the embankment but a significant anomaly appeared at its base. Drilling revealed that the anomaly was a permeable layer of coarse-grained gravel. Further study of the anomaly identified this area of the embankment as being at risk

**4. Georadar for monitoring soil contamination (Sylwia Tomecka-Suchoń)** 

Georadar methods (GPR) are a valuable tool in detection and remediation of soil pollution due to their sensitivity to the electrical properties of subsurface materials. Hydrocarbons and salts leaking from containment at industrial sites are the two most common classes of soil pollutants. Hydrocarbon contamination is usually a consequence of leakage from collecting tanks and pipelines used in various petroleum operations. Mining operations with their sizeable holding ponds, debris, and chemical dumps can also pollute soil and ground water. Identification and remediation of hazardous chemicals is complicated by the fact that waste material can be distributed at various depths, and may migrate in any direction, at any speed. The GPR method is a particularly sensitive method under conditions of pore water saturation, due to the high relative dielectric permittivity of water. In saturated soils for example, GPR can be used to perform real time monitoring of contaminant plume migration. As with applications described in the previous section, GPR is non-invasive and allows parties involved to forgo the more expensive activities of digging wells and trenches (Gołębiowski et al., 2010a). Two case studies of GPR use in soil contamination studies are described below. These include both types of liquid pollutants: low-conductivity hydrocarbons (a fuel station in Krakow) and high-conductivity chemical solutions (former

The conductivity of soil media depends on 1) the relative volumes of the components (e.g., articles, water, air) and 2) their conductivity. The former volumetric factors will in turn affect the soil's overall dielectric permittivity. Owing to its very high dielectric permittivity, (εr =81) water content greatly facilitates imaging of Light Non-Aqueous Phase Liquid (LNAPL) contamination in sandy subsurface materials. Hydrocarbons have a low εr (Table 2.1). Detection of a LNAPL plume in dry, aerated soil horizons is difficult, but in water saturated soils, phase separation occurs between the LNAPL plume and uncontaminated surroundings. Differences in permittivity mean that higher concentrations of immiscible hydrocarbons (e.g. gasoline) can enhance detection and mapping of plumes. LNAPL plumes tend migrate in a downward direction due to gravity and capillary action. When a plume reaches the capillary fringe zone, contrasting electromagnetic properties between

contaminants and the surrounding medium resolve the plume in GPR echograms.

A GPR survey was carried out at a fuel station in Krakow (Gołębiowski et al., 2010a). The study area was contaminated by repeated spills where hydrocarbon mixtures soaked and infiltrated the ground surface. GPR measurements were carried out in a constant-offset

of hydraulic puncture during high water events.

military installation at B-S (Poland), and D waste dump).

**4.1 Low-conductivity hydrocarbon contamination** 

Fig. 3.3 Echogram in the Saint Peter and Saint Paul church in Krakow. RAMAC/CU II apparatus, shielded 250 MHz antennas.

GPR is highly effective in locating zones of weakness in river embankments. This method is useful for performing qualitative investigations of the embankment's inner structure and basal materials. GPR measurement profiles are best carried out along the rim of an embankment, along inner and outer shelves, and also within intervening areas between the river and the embankment. GPR signals can render the subsurface in sufficiently high resolution as to show relatively small anomalies in the embankment's subsurface structure. Varying degrees of cementation for example can be detected by GPR. Lower frequency antennas are generally used to query an embankment's basal structure.

Fig. 3.4 Echogram of the Odra river embankment. RAMAC/GPR apparatus, unshielded 200 MHz antennas.

Frequently, river embankments are composed of impermeable clays and silts, materials that have relatively low resistivity. These materials attenuate EM signals, allowing GPR surveys to penetrate only the shallowest zones of the subsurface. GPR signals can also be obscured by high water, flooding, or heavy rains which saturate embankment material, and thus attenuate EM waves (Nguyen et al., 2005). Trees growing on or near the embankment can also interfere with signals, obscuring recognition of heterogeneities.

In spite of these limitations, the high measurement efficiency and low cost of GPR relative to other methods recommend it as a first order survey strategy for assessing embankments. GPR anomalies can then be confirmed with second order geotechnical methods. Used in this

Fig. 3.3 Echogram in the Saint Peter and Saint Paul church in Krakow. RAMAC/CU II

antennas are generally used to query an embankment's basal structure.

also interfere with signals, obscuring recognition of heterogeneities.

GPR is highly effective in locating zones of weakness in river embankments. This method is useful for performing qualitative investigations of the embankment's inner structure and basal materials. GPR measurement profiles are best carried out along the rim of an embankment, along inner and outer shelves, and also within intervening areas between the river and the embankment. GPR signals can render the subsurface in sufficiently high resolution as to show relatively small anomalies in the embankment's subsurface structure. Varying degrees of cementation for example can be detected by GPR. Lower frequency

Fig. 3.4 Echogram of the Odra river embankment. RAMAC/GPR apparatus, unshielded

Frequently, river embankments are composed of impermeable clays and silts, materials that have relatively low resistivity. These materials attenuate EM signals, allowing GPR surveys to penetrate only the shallowest zones of the subsurface. GPR signals can also be obscured by high water, flooding, or heavy rains which saturate embankment material, and thus attenuate EM waves (Nguyen et al., 2005). Trees growing on or near the embankment can

In spite of these limitations, the high measurement efficiency and low cost of GPR relative to other methods recommend it as a first order survey strategy for assessing embankments. GPR anomalies can then be confirmed with second order geotechnical methods. Used in this

apparatus, shielded 250 MHz antennas.

200 MHz antennas.

way, GPR can optimize the more expensive assessment techniques, limiting their use to zones where GPR anomalies are apparent. The echogram shown in Fig. 3.4 was recorded on the Odra river embankment using 200 MHz unshielded antennas. The height of the embankment along the section analyzed was about 4 m. The goal of the investigation was to identify heterogeneities within the embankment and its basal materials. No meaningful heterogeneities were observed within the embankment but a significant anomaly appeared at its base. Drilling revealed that the anomaly was a permeable layer of coarse-grained gravel. Further study of the anomaly identified this area of the embankment as being at risk of hydraulic puncture during high water events.

## **4. Georadar for monitoring soil contamination (Sylwia Tomecka-Suchoń)**

Georadar methods (GPR) are a valuable tool in detection and remediation of soil pollution due to their sensitivity to the electrical properties of subsurface materials. Hydrocarbons and salts leaking from containment at industrial sites are the two most common classes of soil pollutants. Hydrocarbon contamination is usually a consequence of leakage from collecting tanks and pipelines used in various petroleum operations. Mining operations with their sizeable holding ponds, debris, and chemical dumps can also pollute soil and ground water. Identification and remediation of hazardous chemicals is complicated by the fact that waste material can be distributed at various depths, and may migrate in any direction, at any speed. The GPR method is a particularly sensitive method under conditions of pore water saturation, due to the high relative dielectric permittivity of water. In saturated soils for example, GPR can be used to perform real time monitoring of contaminant plume migration. As with applications described in the previous section, GPR is non-invasive and allows parties involved to forgo the more expensive activities of digging wells and trenches (Gołębiowski et al., 2010a). Two case studies of GPR use in soil contamination studies are described below. These include both types of liquid pollutants: low-conductivity hydrocarbons (a fuel station in Krakow) and high-conductivity chemical solutions (former military installation at B-S (Poland), and D waste dump).

## **4.1 Low-conductivity hydrocarbon contamination**

The conductivity of soil media depends on 1) the relative volumes of the components (e.g., articles, water, air) and 2) their conductivity. The former volumetric factors will in turn affect the soil's overall dielectric permittivity. Owing to its very high dielectric permittivity, (εr =81) water content greatly facilitates imaging of Light Non-Aqueous Phase Liquid (LNAPL) contamination in sandy subsurface materials. Hydrocarbons have a low εr (Table 2.1). Detection of a LNAPL plume in dry, aerated soil horizons is difficult, but in water saturated soils, phase separation occurs between the LNAPL plume and uncontaminated surroundings. Differences in permittivity mean that higher concentrations of immiscible hydrocarbons (e.g. gasoline) can enhance detection and mapping of plumes. LNAPL plumes tend migrate in a downward direction due to gravity and capillary action. When a plume reaches the capillary fringe zone, contrasting electromagnetic properties between contaminants and the surrounding medium resolve the plume in GPR echograms.

A GPR survey was carried out at a fuel station in Krakow (Gołębiowski et al., 2010a). The study area was contaminated by repeated spills where hydrocarbon mixtures soaked and infiltrated the ground surface. GPR measurements were carried out in a constant-offset

Geophysics in Near-Surface Investigations 59

contaminated versus uncontaminated areas can be differentiated according to signal amplitudes in echograms (Marcak & Tomecka-Suchoń, 2010). The higher conductivity of some pollutants leads to greater attenuation and in echograms the absence of reflexes may be expected. Real time monitoring at the B-S military site revealed immediate changes in the GPR signals following injection of a high-conductivity NaCl solution into a zone of concentrated hydrocarbon contamination. The salt solution (5 kg of NaCl in 30 liters of water) was injected into a well located at point C along the profile (Fig. 4.3). Twenty liters of gasoline were simultaneously poured into the ground surface at point B. The injection well passes through sands of mixed grain size, overlying a clay layer. The clay layer is impermeable, separating the contamination experiment from the top of the water table

GPR recordings were performed within the aeration zone of a lossless geological medium using high resolution antennas (500 and 800 MHz) at low depth, and operating in timemonitoring mode in order to track migration of contaminants (Fig. 4.3). A mean velocity of Vmean= 10 cm/ns, estimated from well core materials (i.e. dry sands) was used for timedepth conversions. GPR results show that attenuation increases dramatically around the borehole such that the lateral boundary of the contaminant plume was easily identified. The depth of contamination however was more difficult to locate due to the high conductivity and high attenuation of the NaCl solution. GPR time-monitoring of the gasoline spill showed that the horizontal length of the salt-contaminated area (~30 - 32.5 m; Fig. 4.3a-d) exceeded the original area of the spill, and did not change significantly. Twenty four hours later the contaminant plume had moved down on the profile and a hyperbola located at 30.5

m which probably originated from root reappeared on GPR echograms (Fig. 4.3d).

Different areas of the same water table can be contaminated from different local point sources, leading to varying levels of contamination along a given horizontal transect. In these situations, a single attenuation coefficient cannot be used to map the contamination. Changes detected in the coefficients themselves however may be used to locate boundaries of various contaminated areas. For porous media, reflection coefficients for EM waves depend not only on pore fluid type, but also on its degree of saturation (Schön, 1996). The effects of conductivity variation (ionic content of pore fluids) are evident in the GPR profile of an area of D waste dump contaminated by a high-conductivity solution, shown in Fig. 4.4. GPR data was winnowed according to an echogram time window of 20 and 35 ns to optimize recognition of the near surface aquifer and the horizon of salt

The energy component of the GPR signal reflected from the salt-contaminated horizon was subjected to a Hilbert transform, and then normalized along both horizontal and vertical axes to determine the dominant trend of energy changes (Fig. 4.4a, blue dots). Energy changes were then projected, using a cubic polynomial function, along the entire profile (Fig. 4.4a, red line). GPR results were merged with hydrological data on the basis of pore fluid conductivity (green line). Near the piezometer Pz-2, concentrated pore fluids corresponded to high energy GPR signals. To the right of Pz-2, GPR signal energy decreases rapidly until well site S-640A. The trend in energy change (Fig. 4.4a, red line) corresponds to

trends in water conductivity change (Fig. 4.4a, green line).

at ~9 m depth.

contamination.

reflection mode along two lines that transected the site. A RAMAC/GPR device was used with 200 MHz antennas. Traces were collected every 0.05 m to generate a total of 64 stacked signals. All echograms were processed with the ReflexW program using phase correlation, time zero correction, amplitude declipping, dewowing, DC-shift, background removal, gain, the Butterworth filter, and smoothing procedures. Identification of contaminated zones was based on energy distribution analysis, counted from Hilbert transformation of raw signals. Areas where pockets of high energy (Fig. 4.1) overlapped anomalies in power spectra were interpreted as contaminated zones (Fig. 4.2). The energy envelope was normalized to the maximum value of the direct air wave.

Analysis of well cores (Fig. 4.1) found neither anisotropy nor other compositional anomalies in the subsurface of the study area. Parts of high energy in the echograms (violet area between 1 m and 2.5 m in Fig. 4.1) were therefore interpreted as contaminated zones. Analysis of GPR power spectra from the profile (Fig. 4.2) showed phase shifts and other irregularities that likewise indicated high concentrations of hydrocarbons. Power spectra irregularities correlated well with the zones of high energy in areas located at 15 and 32 m along the transect (Figs 4.1 and 4.2). The location of a high energy zone at ~2.5 m depth highlighted depression or downward bowing of the water table due to the gravitational pressure of the hydrocarbon plume.

Fig. 4.1 GPR energy distribution (Gołębiowski et al., 2010a)

Fig. 4.2 GPR power spectra (Gołębiowski et al., 2010a)

#### **4.2 High-conductivity chemical contamination**

GPR can be used to detect subsurface chemical contamination through analysis of electrical conductivity profiles. Conductivity strongly affects attenuation of EM waves such that

reflection mode along two lines that transected the site. A RAMAC/GPR device was used with 200 MHz antennas. Traces were collected every 0.05 m to generate a total of 64 stacked signals. All echograms were processed with the ReflexW program using phase correlation, time zero correction, amplitude declipping, dewowing, DC-shift, background removal, gain, the Butterworth filter, and smoothing procedures. Identification of contaminated zones was based on energy distribution analysis, counted from Hilbert transformation of raw signals. Areas where pockets of high energy (Fig. 4.1) overlapped anomalies in power spectra were interpreted as contaminated zones (Fig. 4.2). The energy envelope was normalized to the

Analysis of well cores (Fig. 4.1) found neither anisotropy nor other compositional anomalies in the subsurface of the study area. Parts of high energy in the echograms (violet area between 1 m and 2.5 m in Fig. 4.1) were therefore interpreted as contaminated zones. Analysis of GPR power spectra from the profile (Fig. 4.2) showed phase shifts and other irregularities that likewise indicated high concentrations of hydrocarbons. Power spectra irregularities correlated well with the zones of high energy in areas located at 15 and 32 m along the transect (Figs 4.1 and 4.2). The location of a high energy zone at ~2.5 m depth highlighted depression or downward bowing of the water table due to the gravitational

GPR can be used to detect subsurface chemical contamination through analysis of electrical conductivity profiles. Conductivity strongly affects attenuation of EM waves such that

maximum value of the direct air wave.

pressure of the hydrocarbon plume.

Fig. 4.1 GPR energy distribution (Gołębiowski et al., 2010a)

Fig. 4.2 GPR power spectra (Gołębiowski et al., 2010a)

**4.2 High-conductivity chemical contamination** 

contaminated versus uncontaminated areas can be differentiated according to signal amplitudes in echograms (Marcak & Tomecka-Suchoń, 2010). The higher conductivity of some pollutants leads to greater attenuation and in echograms the absence of reflexes may be expected. Real time monitoring at the B-S military site revealed immediate changes in the GPR signals following injection of a high-conductivity NaCl solution into a zone of concentrated hydrocarbon contamination. The salt solution (5 kg of NaCl in 30 liters of water) was injected into a well located at point C along the profile (Fig. 4.3). Twenty liters of gasoline were simultaneously poured into the ground surface at point B. The injection well passes through sands of mixed grain size, overlying a clay layer. The clay layer is impermeable, separating the contamination experiment from the top of the water table at ~9 m depth.

GPR recordings were performed within the aeration zone of a lossless geological medium using high resolution antennas (500 and 800 MHz) at low depth, and operating in timemonitoring mode in order to track migration of contaminants (Fig. 4.3). A mean velocity of Vmean= 10 cm/ns, estimated from well core materials (i.e. dry sands) was used for timedepth conversions. GPR results show that attenuation increases dramatically around the borehole such that the lateral boundary of the contaminant plume was easily identified. The depth of contamination however was more difficult to locate due to the high conductivity and high attenuation of the NaCl solution. GPR time-monitoring of the gasoline spill showed that the horizontal length of the salt-contaminated area (~30 - 32.5 m; Fig. 4.3a-d) exceeded the original area of the spill, and did not change significantly. Twenty four hours later the contaminant plume had moved down on the profile and a hyperbola located at 30.5 m which probably originated from root reappeared on GPR echograms (Fig. 4.3d).

Different areas of the same water table can be contaminated from different local point sources, leading to varying levels of contamination along a given horizontal transect. In these situations, a single attenuation coefficient cannot be used to map the contamination. Changes detected in the coefficients themselves however may be used to locate boundaries of various contaminated areas. For porous media, reflection coefficients for EM waves depend not only on pore fluid type, but also on its degree of saturation (Schön, 1996). The effects of conductivity variation (ionic content of pore fluids) are evident in the GPR profile of an area of D waste dump contaminated by a high-conductivity solution, shown in Fig. 4.4. GPR data was winnowed according to an echogram time window of 20 and 35 ns to optimize recognition of the near surface aquifer and the horizon of salt contamination.

The energy component of the GPR signal reflected from the salt-contaminated horizon was subjected to a Hilbert transform, and then normalized along both horizontal and vertical axes to determine the dominant trend of energy changes (Fig. 4.4a, blue dots). Energy changes were then projected, using a cubic polynomial function, along the entire profile (Fig. 4.4a, red line). GPR results were merged with hydrological data on the basis of pore fluid conductivity (green line). Near the piezometer Pz-2, concentrated pore fluids corresponded to high energy GPR signals. To the right of Pz-2, GPR signal energy decreases rapidly until well site S-640A. The trend in energy change (Fig. 4.4a, red line) corresponds to trends in water conductivity change (Fig. 4.4a, green line).

Geophysics in Near-Surface Investigations 61

Fig. 4.4 Results of GPR investigations of an area surrounding a mining dump; a) energy distributions of GPR signals along the measurement profile and changes of ground water

perpendicular stripes are zones of high interference that were removed from GPR signals

Micro-scale gravity measurements are the newest type of the gravimetric investigations in geophysics. Anomalies measured in the microgravity range are much smaller than those measured by traditional gravimetric methods. Microgravity amplitudes are only few times higher than average measurement errors for typical gravimetric devices. The newest generation of gravity meters have measurement accuracies of ~5 x 10-8 m/s2 (0.005 mGal). Gravity field values for a given object depend primarily on how its bulk density is spatially distributed. For subsurface materials, shapes, sizes and the depth of objects may effect gravity measurements, but other petrophysical properties do not. Microgravity measurements can be influenced by strong winds or ground vibrations from traffic, both of which can lead to large errors. Special methods can be used to reduce the microgravity measurement error and enhance a signal/error ratio. These include shortening the time intervals between measurements to eliminate instrument drift (i.e. measurement intervals of no more than 1-3 hrs), performing replicate measurements to reduce random measurement error, and using relatively short distances between measurement stations in order to image anomalies as precisely as possible. Gravimetric modelling of the geological or human - made

conductivity according to hydrogeological results; b) modified echogram - white

**5. Microgravity investigations (Sławomir Porzucek)** 

(Gołębiowski et al., 2010b)

Fig. 4.3 Echograms (800 MHz antennas) from the B-S site a) before simulated gasoline spill, b) immediately after the spill, c) 3 hours after d), 24 hours after (Gołębiowski *et al.*, 2010a)

Fig. 4.3 Echograms (800 MHz antennas) from the B-S site a) before simulated gasoline spill, b) immediately after the spill, c) 3 hours after d), 24 hours after (Gołębiowski *et al.*, 2010a)

Fig. 4.4 Results of GPR investigations of an area surrounding a mining dump; a) energy distributions of GPR signals along the measurement profile and changes of ground water conductivity according to hydrogeological results; b) modified echogram - white perpendicular stripes are zones of high interference that were removed from GPR signals (Gołębiowski et al., 2010b)

## **5. Microgravity investigations (Sławomir Porzucek)**

Micro-scale gravity measurements are the newest type of the gravimetric investigations in geophysics. Anomalies measured in the microgravity range are much smaller than those measured by traditional gravimetric methods. Microgravity amplitudes are only few times higher than average measurement errors for typical gravimetric devices. The newest generation of gravity meters have measurement accuracies of ~5 x 10-8 m/s2 (0.005 mGal). Gravity field values for a given object depend primarily on how its bulk density is spatially distributed. For subsurface materials, shapes, sizes and the depth of objects may effect gravity measurements, but other petrophysical properties do not. Microgravity measurements can be influenced by strong winds or ground vibrations from traffic, both of which can lead to large errors. Special methods can be used to reduce the microgravity measurement error and enhance a signal/error ratio. These include shortening the time intervals between measurements to eliminate instrument drift (i.e. measurement intervals of no more than 1-3 hrs), performing replicate measurements to reduce random measurement error, and using relatively short distances between measurement stations in order to image anomalies as precisely as possible. Gravimetric modelling of the geological or human - made

Geophysics in Near-Surface Investigations 63

m, and occurring at depths of several meters up to several tens of meters. Archival maps provided the approximate location of these shafts. The survey was conducted using a main grid of 50 x 50 m, with stations spaced every 5 m. Additional stations were implemented within the main grid, creating a secondary 17.5 x 17.5 m grid with stations every 2.5 m.

Fig. 5.1 Microgravity residual anomalies in the eastern part of the Wieliczka salt mine

Fig. 5.2 Bouguer anomaly in the area where the shaft was suspected to be located; irregular positioning of stations is due to terrain, as well as surface and underground infrastructure

structures under investigation can assist in selecting the optimal distance interval between measurement stations. The height of the tripod on which the instrument is mounted should be measured as precisely as possible. An error of 0.03 m in tripod height estimation can change gravity measurements by as much as 0.01 mGal. Station location and terrain corrections are also extremely important. Improper processing of microgravity data (e.g. poorly chosen filters and so forth) may mask subsurface anomalies or create artefacts in the data. A poorly executed microgravity survey may reveal little more than the obvious Bouguer anomalies. For a carefully executed survey, microgravity anomalies with amplitudes of as little as 1.5 times the measurement error can be identified from as few as four stations, with a 99.9% confidence level (Liu, 2007). Microgravity methods are highly sensitive to the density of subsurface objects and thus are generally used for investigating small objects at relatively shallow depths. Microgravimetric information combined with borehole data can provide 2- and 3-D density models of the subsurface. Bulk density differences apparent in microgravity data can be to used locate and identify loosened debris, such as slag heaps, salts or other industrial waste materials. Microgravity surveys are also useful in searching for natural cavities caused by karst processes. In both active and abandoned subsurface mine sites, microgravity surveys can show underground shafts, chambers and other low density features. Microgravity surveys can thus identify risks to surface activities posed by old, closed underground workings and their potential deformation.

A microgravity survey of the Wieliczka salt mine provides an example of the application described above. A 600 x 350 m area was surveyed for microgravity anomalies to indentify zones in the rock mass that posed a structural risk to surface activities (Madej *et al*., 2001). Residual anomalies (Fig. 5.1) reflected both the geological structure of the subsurface and structures related to former mining operations. Negative concentric anomalies revealed chambers excavated by historical mining operations. The negative, concentric anomaly denoted as '1' in Fig. 5.1, was the result of weakened materials around an abandoned 16th century structure, the Lois shaft. Negative anomalies 2, 3 and 4 (Fig. 5.1) were interpreted as the main excavation area of the 16th century operation. Collapsing mine chambers and corresponding surface subsidence were recorded as early as the 16th century, indicating an on-going risk to modern surface activities. Microgravimetric surveys and modelling based on historical information enabled location of the weakened zones between level I of the mine and the surface. The negative, concentric anomalies 5 and 6 (Fig. 5.1) corresponded to excavation chambers located on level II. The relatively low amplitudes of these anomalies indicated insignificant dewatering of the strata above the chambers, and thus a slower than expected advance of the weakened zones towards the surface. A WWN-EES trending negative anomaly 7; (Fig. 5.1) is situated near a local road. This anomaly corresponded to a system of chambers that was undergoing uplift and deformation. The collapse of roof material caused structurally weakened zones above the chambers and posed significant risk to the road.

Microgravity surveys have also been conducted near a calamine deposit that was exploited during the first half of the 19th century. These investigations sought to identify shafts that posed risk to planned construction of a nearby road. The calamine deposit occurred in a Lower Triassic dolomite unit directly overlain by relatively thin (< 1 m) Quaternary deposits. The exploited calamine horizon was about 2 m thick, and located at depths of ~20- 30 m. The mining operation consisted of a system of narrow shafts having diameters of 1.5

structures under investigation can assist in selecting the optimal distance interval between measurement stations. The height of the tripod on which the instrument is mounted should be measured as precisely as possible. An error of 0.03 m in tripod height estimation can change gravity measurements by as much as 0.01 mGal. Station location and terrain corrections are also extremely important. Improper processing of microgravity data (e.g. poorly chosen filters and so forth) may mask subsurface anomalies or create artefacts in the data. A poorly executed microgravity survey may reveal little more than the obvious Bouguer anomalies. For a carefully executed survey, microgravity anomalies with amplitudes of as little as 1.5 times the measurement error can be identified from as few as four stations, with a 99.9% confidence level (Liu, 2007). Microgravity methods are highly sensitive to the density of subsurface objects and thus are generally used for investigating small objects at relatively shallow depths. Microgravimetric information combined with borehole data can provide 2- and 3-D density models of the subsurface. Bulk density differences apparent in microgravity data can be to used locate and identify loosened debris, such as slag heaps, salts or other industrial waste materials. Microgravity surveys are also useful in searching for natural cavities caused by karst processes. In both active and abandoned subsurface mine sites, microgravity surveys can show underground shafts, chambers and other low density features. Microgravity surveys can thus identify risks to surface activities posed by old, closed underground workings and their potential

A microgravity survey of the Wieliczka salt mine provides an example of the application described above. A 600 x 350 m area was surveyed for microgravity anomalies to indentify zones in the rock mass that posed a structural risk to surface activities (Madej *et al*., 2001). Residual anomalies (Fig. 5.1) reflected both the geological structure of the subsurface and structures related to former mining operations. Negative concentric anomalies revealed chambers excavated by historical mining operations. The negative, concentric anomaly denoted as '1' in Fig. 5.1, was the result of weakened materials around an abandoned 16th century structure, the Lois shaft. Negative anomalies 2, 3 and 4 (Fig. 5.1) were interpreted as the main excavation area of the 16th century operation. Collapsing mine chambers and corresponding surface subsidence were recorded as early as the 16th century, indicating an on-going risk to modern surface activities. Microgravimetric surveys and modelling based on historical information enabled location of the weakened zones between level I of the mine and the surface. The negative, concentric anomalies 5 and 6 (Fig. 5.1) corresponded to excavation chambers located on level II. The relatively low amplitudes of these anomalies indicated insignificant dewatering of the strata above the chambers, and thus a slower than expected advance of the weakened zones towards the surface. A WWN-EES trending negative anomaly 7; (Fig. 5.1) is situated near a local road. This anomaly corresponded to a system of chambers that was undergoing uplift and deformation. The collapse of roof material caused structurally weakened zones above the chambers and posed significant risk

Microgravity surveys have also been conducted near a calamine deposit that was exploited during the first half of the 19th century. These investigations sought to identify shafts that posed risk to planned construction of a nearby road. The calamine deposit occurred in a Lower Triassic dolomite unit directly overlain by relatively thin (< 1 m) Quaternary deposits. The exploited calamine horizon was about 2 m thick, and located at depths of ~20- 30 m. The mining operation consisted of a system of narrow shafts having diameters of 1.5

deformation.

to the road.

m, and occurring at depths of several meters up to several tens of meters. Archival maps provided the approximate location of these shafts. The survey was conducted using a main grid of 50 x 50 m, with stations spaced every 5 m. Additional stations were implemented within the main grid, creating a secondary 17.5 x 17.5 m grid with stations every 2.5 m.

Fig. 5.1 Microgravity residual anomalies in the eastern part of the Wieliczka salt mine

Fig. 5.2 Bouguer anomaly in the area where the shaft was suspected to be located; irregular positioning of stations is due to terrain, as well as surface and underground infrastructure

Geophysics in Near-Surface Investigations 65

Fig. 5.4 a) Bouguer anomaly A, b) result of gravimetric modelling of the A anomaly

Fig. 5.5 a) Bouguer anomaly B, b) result of gravimetric modelling of the B anomaly

Microgravity can also be applied to investigate embankments and earthen dams that are at risk of water leakage and possible failure. Figures 5.4 and 5.5 show results from microgravity surveys conducted to detect zones of weakness along embankments that surround an underground water reservoir. Surveys were carried out along the embankment using a measurement interval of 2 m. A terrain correction was applied to the data to address irregularities in the topography of the study area. The survey revealed some relatively negative anomalies, and identified several zones of weakness within the embankment. The anomaly labeled A is based on 14 measurement stations and had an amplitude of about 0.03 mGal (Fig. 5.4a), a factor of 3 times higher than the overall measurement error of the study. The next anomaly was a separate residual anomaly, which was then modelled (results shown in Fig. 5.4b). One zone identified in the survey was found to have a density of 0.30 Mg/m3 less than the average density of surrounding material in the embankment. This zone begins 2 m below the surface, has a thickness of about 3 m, and contains piping that carries water from the reservoir. Poor sealing of the pipes had likely caused leakage and soil

The number of stations used in this two grid survey strategy provided the subsurface resolution necessary to locate several shafts and surrounding zones of weakness that were at risk of collapse. The highlighted area shown along with the Bouguer anomalies in Fig. 5.2, is a 40 x 20 m zone that exhibited a distinct, negative gravity anomaly. This area encircled a smaller anomaly with a diameter of about 10 m. Old mining maps indicated uncertainty concerning the location of a shaft that might have been moved from its original position in the central part of the study area. The microgravimetric anomalies clearly identified the shaft as the smaller, encircled anomaly. The larger oval-shaped anomaly surrounding the shaft was interpreted as a weakened zone, and identified as a possible risk to road construction and other surface activities.

Fig. 5.3 Gravity anomalies and corrections over the subsurface drift

Microgravity surveys conducted prior to land development projects have also confirmed usefulness of the technique in risk mitigation. Results shown in figure 5.3 were collected from an area overlying a sealed mine drift, tunneled out of the local sandstone formation at the end of the 19th century. The tunnel had a total length of about 15.5 m, a triangular crosssection with a 3-5 m base, and a height of 3-4 m. The microgravity measurement stations were positioned across the drift at intervals of 1-4 m. The drift is visible as a relatively negative anomaly (Fig. 5.3, ΔgB) but the amplitude of the anomaly is difficult to determine due to the bulk density properties of the surrounding media. Excavation of the entrance to the drift increased the difficulty in detecting the structure by gravimetric methods. To minimize the impact of irregularities in the overlying terrain, the triangulation method was used to calculate a terrain correction, gt (Wójcicki, 1993). Investigators also corrected for three dimensional gravitational effects of the drift, gd, using geodesy measurements to subtract the drift from the local gravity anomaly. Together, these corrections yielded a final anomaly, gBdt (Fig. 5.3). The final anomaly was based on points 8–9, and bore some similarity to the initial distribution. The shape of the anomaly however did not correspond to the shape of the applied corrections. The final anomaly was thus interpreted as a fractured zone, weakened by earlier construction of the drift. These findings demonstrated that the subsurface was not entirely stable.

The number of stations used in this two grid survey strategy provided the subsurface resolution necessary to locate several shafts and surrounding zones of weakness that were at risk of collapse. The highlighted area shown along with the Bouguer anomalies in Fig. 5.2, is a 40 x 20 m zone that exhibited a distinct, negative gravity anomaly. This area encircled a smaller anomaly with a diameter of about 10 m. Old mining maps indicated uncertainty concerning the location of a shaft that might have been moved from its original position in the central part of the study area. The microgravimetric anomalies clearly identified the shaft as the smaller, encircled anomaly. The larger oval-shaped anomaly surrounding the shaft was interpreted as a weakened zone, and identified as a possible risk to road

construction and other surface activities.

that the subsurface was not entirely stable.

Fig. 5.3 Gravity anomalies and corrections over the subsurface drift

Microgravity surveys conducted prior to land development projects have also confirmed usefulness of the technique in risk mitigation. Results shown in figure 5.3 were collected from an area overlying a sealed mine drift, tunneled out of the local sandstone formation at the end of the 19th century. The tunnel had a total length of about 15.5 m, a triangular crosssection with a 3-5 m base, and a height of 3-4 m. The microgravity measurement stations were positioned across the drift at intervals of 1-4 m. The drift is visible as a relatively negative anomaly (Fig. 5.3, ΔgB) but the amplitude of the anomaly is difficult to determine due to the bulk density properties of the surrounding media. Excavation of the entrance to the drift increased the difficulty in detecting the structure by gravimetric methods. To minimize the impact of irregularities in the overlying terrain, the triangulation method was used to calculate a terrain correction, gt (Wójcicki, 1993). Investigators also corrected for three dimensional gravitational effects of the drift, gd, using geodesy measurements to subtract the drift from the local gravity anomaly. Together, these corrections yielded a final anomaly, gBdt (Fig. 5.3). The final anomaly was based on points 8–9, and bore some similarity to the initial distribution. The shape of the anomaly however did not correspond to the shape of the applied corrections. The final anomaly was thus interpreted as a fractured zone, weakened by earlier construction of the drift. These findings demonstrated

Fig. 5.4 a) Bouguer anomaly A, b) result of gravimetric modelling of the A anomaly

Fig. 5.5 a) Bouguer anomaly B, b) result of gravimetric modelling of the B anomaly

Microgravity can also be applied to investigate embankments and earthen dams that are at risk of water leakage and possible failure. Figures 5.4 and 5.5 show results from microgravity surveys conducted to detect zones of weakness along embankments that surround an underground water reservoir. Surveys were carried out along the embankment using a measurement interval of 2 m. A terrain correction was applied to the data to address irregularities in the topography of the study area. The survey revealed some relatively negative anomalies, and identified several zones of weakness within the embankment. The anomaly labeled A is based on 14 measurement stations and had an amplitude of about 0.03 mGal (Fig. 5.4a), a factor of 3 times higher than the overall measurement error of the study. The next anomaly was a separate residual anomaly, which was then modelled (results shown in Fig. 5.4b). One zone identified in the survey was found to have a density of 0.30 Mg/m3 less than the average density of surrounding material in the embankment. This zone begins 2 m below the surface, has a thickness of about 3 m, and contains piping that carries water from the reservoir. Poor sealing of the pipes had likely caused leakage and soil

Geophysics in Near-Surface Investigations 67

concentrate from soil samples revealed high concentrations of magnetite. Maghemite and iron sulfides were also identified. Subsequent studies have revealed a significant increase in the magnetic susceptibility of soil closest to the steel plant in the last four years (Rosowiecka & Nawrocki, 2010). Land near the steel plant is composed of fertile soil (chernozem, brown soils, loess) and is under cultivation. Magnetometry is especially suitable for further

Fig. 6.1 Map of the magnetic susceptibility of topsoil in the vicinity of a Krakow steel plant

A magnetic susceptibility study of historical sequence layers underlying the Main Market Square in Krakow revealed information about past human activity in the area. Subsurface samples collected during an archaeological excavation were analyzed in the laboratory using a MS2 meter and MS2B sensor. Analysis revealed contrasts in the magnetic susceptibility of underlying sands and overlying layers affected by anthropogenic activities (Fig. 6.2). Anthropogenic layers, which contained dark particles, fragments of bricks and sediment filled void space, often showed greater magnetic susceptibility. Overall, the

Fig. 6.2 Magnetic susceptibility of samples (χ · 10-8 m3kg-1) of historical sequence layers from

**6.3 Magnetic susceptibility of historical sequence layers** 

historical layers had relatively weak magnetic properties.

the Main Market Square in Krakow.

effective monitoring of the surrounding area.

suffusion, resulting in this lower density zone. The second anomaly marked B had a much smaller amplitude (Fig 5.5a) and was located beneath the causeway crossing the reservoir. The anomaly is based on 6 - 7 measurement stations, and was identified with a 99.99% confidence level (Liu, 2007).

The B anomaly was modelled (Fig. 5.5b) under the assumptions that it was a zone of poorly consolidated material that extended from the bottom of the reservoir to a depth of about 2 m below the surface. The top of the zone of weakness was also analyzed using GPR echograms, allowing for a consistent and unambiguous model solution. The left edge of the unconsolidated zone corresponded to the location of an HDPE reservoir liner, embedded within the slope of the main embankment. In this case, microgravimetric identification of the zone of weakness demonstrated poor construction of the embankment in the study area.

## **6. Magnetometry in the investigation of surface materials (Anna Wojas)**

## **6.1 Origin of magnetic particles**

In near surface investigations of the Earth's crust (shallow rock formations and soils), magnetometry can be used to identify magnetic minerals, determine their origin and reveal certain physicochemical processes that may be occurring in the soil. Soil magnetometry measures a given soil's magnetic susceptibility. The technique is commonly used to assess topsoil contamination by heavy metals and to track sources of magnetic particles in urban and industrial areas. Fly ash from incinerators is among others the source of magnetic particles in urban and industrial soils. The magnetic particles derive mostly from combustion of coals containing iron sulfides. High temperature oxidation transforms iron sulfides into magnetite, maghemite, and other less common types of ferrites. Both human activities and natural processes can influence the magnetic susceptibility of soils. Magnetic anomalies can thus be caused by industrial operations, weathering of ferrimagnetic and antiferromagnetic minerals, or other pedogenic processes.

#### **6.2 The role of magnetic susceptibility in assessing topsoil contamination**

Surface surveys of the magnetic susceptibility of topsoil were carried out in the vicinity of a major steel operation in Krakow. The soil surrounding the steel plant is contaminated by high levels of Zn, Pb, Cu, Cd and Hg.

Magnetic susceptibility surveys were performed only in the northeastern corner of the steel plant due to limitations in site accessibility. Prevailing winds blow emissions from the plant in a southeasterly direction. A nearby railroad track was also considered as a possible source of magnetic particles in study area soils. *In situ* measurements of magnetic susceptibility were performed using a MS2 meter and a MS2D sensor (Bartington Co.), working from principles of magnetic induction. Measured magnetic susceptibilities of soil were in the range of 20 to 435 x 10-5 SI units. The highest values were found closest to the steel plant (average values of 215 x 10-5 SI units). These results confirmed fall and nearby deposition of the heavy fraction of emission particles in the soil. The outlying regions of the study area had soil magnetic susceptibility ranging from 20 to 122 x 10-5 SI units, with an average value of 64 x 10-5 SI units (Fig. 6.1). An oval shaped magnetic susceptibility anomaly (isoline 150 x 10-5 SI units) having an easterly directed axis, highlights the direction of prevailing winds and magnetic particle transport. Magnetic and mineralogical analyses of the magnetic

suffusion, resulting in this lower density zone. The second anomaly marked B had a much smaller amplitude (Fig 5.5a) and was located beneath the causeway crossing the reservoir. The anomaly is based on 6 - 7 measurement stations, and was identified with a 99.99%

The B anomaly was modelled (Fig. 5.5b) under the assumptions that it was a zone of poorly consolidated material that extended from the bottom of the reservoir to a depth of about 2 m below the surface. The top of the zone of weakness was also analyzed using GPR echograms, allowing for a consistent and unambiguous model solution. The left edge of the unconsolidated zone corresponded to the location of an HDPE reservoir liner, embedded within the slope of the main embankment. In this case, microgravimetric identification of the zone of weakness demonstrated poor construction of the embankment in the study area.

**6. Magnetometry in the investigation of surface materials (Anna Wojas)** 

**6.2 The role of magnetic susceptibility in assessing topsoil contamination** 

Surface surveys of the magnetic susceptibility of topsoil were carried out in the vicinity of a major steel operation in Krakow. The soil surrounding the steel plant is contaminated by

Magnetic susceptibility surveys were performed only in the northeastern corner of the steel plant due to limitations in site accessibility. Prevailing winds blow emissions from the plant in a southeasterly direction. A nearby railroad track was also considered as a possible source of magnetic particles in study area soils. *In situ* measurements of magnetic susceptibility were performed using a MS2 meter and a MS2D sensor (Bartington Co.), working from principles of magnetic induction. Measured magnetic susceptibilities of soil were in the range of 20 to 435 x 10-5 SI units. The highest values were found closest to the steel plant (average values of 215 x 10-5 SI units). These results confirmed fall and nearby deposition of the heavy fraction of emission particles in the soil. The outlying regions of the study area had soil magnetic susceptibility ranging from 20 to 122 x 10-5 SI units, with an average value of 64 x 10-5 SI units (Fig. 6.1). An oval shaped magnetic susceptibility anomaly (isoline 150 x 10-5 SI units) having an easterly directed axis, highlights the direction of prevailing winds and magnetic particle transport. Magnetic and mineralogical analyses of the magnetic

In near surface investigations of the Earth's crust (shallow rock formations and soils), magnetometry can be used to identify magnetic minerals, determine their origin and reveal certain physicochemical processes that may be occurring in the soil. Soil magnetometry measures a given soil's magnetic susceptibility. The technique is commonly used to assess topsoil contamination by heavy metals and to track sources of magnetic particles in urban and industrial areas. Fly ash from incinerators is among others the source of magnetic particles in urban and industrial soils. The magnetic particles derive mostly from combustion of coals containing iron sulfides. High temperature oxidation transforms iron sulfides into magnetite, maghemite, and other less common types of ferrites. Both human activities and natural processes can influence the magnetic susceptibility of soils. Magnetic anomalies can thus be caused by industrial operations, weathering of ferrimagnetic and

confidence level (Liu, 2007).

**6.1 Origin of magnetic particles** 

high levels of Zn, Pb, Cu, Cd and Hg.

antiferromagnetic minerals, or other pedogenic processes.

concentrate from soil samples revealed high concentrations of magnetite. Maghemite and iron sulfides were also identified. Subsequent studies have revealed a significant increase in the magnetic susceptibility of soil closest to the steel plant in the last four years (Rosowiecka & Nawrocki, 2010). Land near the steel plant is composed of fertile soil (chernozem, brown soils, loess) and is under cultivation. Magnetometry is especially suitable for further effective monitoring of the surrounding area.

Fig. 6.1 Map of the magnetic susceptibility of topsoil in the vicinity of a Krakow steel plant

## **6.3 Magnetic susceptibility of historical sequence layers**

A magnetic susceptibility study of historical sequence layers underlying the Main Market Square in Krakow revealed information about past human activity in the area. Subsurface samples collected during an archaeological excavation were analyzed in the laboratory using a MS2 meter and MS2B sensor. Analysis revealed contrasts in the magnetic susceptibility of underlying sands and overlying layers affected by anthropogenic activities (Fig. 6.2). Anthropogenic layers, which contained dark particles, fragments of bricks and sediment filled void space, often showed greater magnetic susceptibility. Overall, the historical layers had relatively weak magnetic properties.

Fig. 6.2 Magnetic susceptibility of samples (χ · 10-8 m3kg-1) of historical sequence layers from the Main Market Square in Krakow.

Geophysics in Near-Surface Investigations 69

small scale faulting and other discontinuities (Goodman & Shi, 1985). Discontinuous deformational features pose the greatest risk to the overlying surface. In solid rocks, seismic surveys can identify fractured zones and potential fault planes as velocity changes in seismic data. A seismic refraction survey of shallow workings at a former mine site within Carboniferous rocks, revealed differing behaviour and degrees of deformation (Table 7.1).

Carboniferous rocks Fractures Hazards to surface

uppermost Carboniferous fractured rocks infiltration hazard

of rocks

rocks

Standard seismic surveys often cannot be conducted at construction sites in developed areas due to their disruptiveness, risks to nearby inhabitants, and accessibility of the area

In these cases, the uppermost body of bedrock beneath the surface can be imaged by seismic refraction tomography (Fig. 7.1). This technique interprets changes in wave velocity during P-wave propagation, in this example, through a near-surface Quaternary deposit (AB and CD in Fig. 7.1; velocity = V0) and underlying bedrock (BC in Fig. 7.1; velocity = V1). P-wave velocity for the Quaternary deposit was much smaller than that of the Carboniferous rocks (V0 < V1), which resulted in a relatively small critical angle of refraction (i < 10o). The near surface Quaternary deposit was only a few meters thick. AB was thus insignificant, and the distance BC was assumed to equal that of AD. The errors for estimating unit thickness and their associated velocity changes were based on critical angle, i, and refractor depth, h.

Given an AD thickness of 100 m, refractor depth h of 2 - 8m, and i equal to 5 - 15o, errors were relatively small, ranging from 0.3 – 4.3 %. In the model described above, the initial arrival time for the refracted signal *TABCD* (traversing ABCD), could be approximated as the time necessary to traverse BC (Dec, 2004). The study used a seismic source located at point B, and receiver at point C, which reduced the initial arrival time to *TAB*, or time necessary to traverse distance AB. Arrival times were used to construct a tomographic image of the Pwave velocity distribution for rocks beneath the Carboniferous surface. Velocity values were

fractured blocks

highly fractured

rocks, fissures

None infiltration hazard at

block contacts

discontinuous deformation hazard

infiltration and discontinuous deformation hazards

high risk of discontinuous deformation

Terrain Category

B >1900

surrounding the site.

Velocity, VP, [m/s] Structure of

block

C 1500-1900 continuous surface of

C1 1500-1900 separate blocks of rocks, displaced

C2 1100-1500 fractured rock, possible

Table 7.1 Seismic categories of Carboniferous rocks

**7.1.1 Seismic surveys in developed areas** 

blocks of rocks

D <1100 destruction of rock mass highly fractured

A >1900 contiguous rock None none

separate blocks of rocks, contiguous rock within a

#### **6.4 Magnetic susceptibility studies of ochra deposit**

Magnetometry studies of an ochre deposit in the Carpathians demonstrate the use of magnetometry in reconnaissance of potential economic deposits (Wojas, 2009). Certain iron oxy- hydroxide minerals weather to an ochre colour, which lends its name to ochre type deposits. These deposits are generally composed of goethite (yellow), hematite (red), and manganese oxides (dark hues). The magnetic susceptibility of the ochre deposit under investigation was attributable to the accumulation of secondary iron sulfides weathered from nearby rocks, and was also influenced by underlying bedrock. Magnetic susceptibility of the area was measured *in situ* using a MS2 meter and MS2F sensor. Low values of magnetic susceptibility (40 to 60 x 10-5 SI units) indicated soil devoid of ochre. The ochre deposit itself was recognizable by its stronger magnetic properties, and wide ranging magnetic susceptibility values. Magnetic susceptibility of yellow and rubiginous ochre ranged from 60 to 200 x 10-5 SI units. Crystalline rocks containing high concentrations of Fe2O3, and appearing as brown ochre had the highest magnetic susceptibility, from 200 to 900 x 10-5 SI units (Kotlarczyk & Ratajczak, 2002) (Fig. 6.3).

Fig. 6.3 Magnetic susceptibility values for an ochre deposit shown with a geological cross section of the deposit (Kotlarczyk & Ratajczak, 2002).

## **7. Environmental seismic investigations (Jerzy Dec)**

#### **7.1 Seismic evaluation of surface risks in developed areas overlying former mine sites**

Closed or abandoned coal mines are often monitored using seismic surveys, which can locate and characterize workings a few meters below the surface, to depths of up to 60 m. Such workings may collapse or destabilize their surroundings, and thus pose a significant risk to surface activities. The stress field surrounding an abandoned working depends on the initial stress affecting the subsurface structure, which is in turn related to the depth and geometry of the working. The stress field can be evident as both continuous and discontinuous deformation (Mange & Kochonov, 1994). Continuous deformation occurs as elastic and plastic strain in layers, which deform but maintain their continuity. Discontinuous deformation manifests as fractures and displacement of bedrock, evident in

Magnetometry studies of an ochre deposit in the Carpathians demonstrate the use of magnetometry in reconnaissance of potential economic deposits (Wojas, 2009). Certain iron oxy- hydroxide minerals weather to an ochre colour, which lends its name to ochre type deposits. These deposits are generally composed of goethite (yellow), hematite (red), and manganese oxides (dark hues). The magnetic susceptibility of the ochre deposit under investigation was attributable to the accumulation of secondary iron sulfides weathered from nearby rocks, and was also influenced by underlying bedrock. Magnetic susceptibility of the area was measured *in situ* using a MS2 meter and MS2F sensor. Low values of magnetic susceptibility (40 to 60 x 10-5 SI units) indicated soil devoid of ochre. The ochre deposit itself was recognizable by its stronger magnetic properties, and wide ranging magnetic susceptibility values. Magnetic susceptibility of yellow and rubiginous ochre ranged from 60 to 200 x 10-5 SI units. Crystalline rocks containing high concentrations of Fe2O3, and appearing as brown ochre had the highest magnetic susceptibility, from 200 to

Fig. 6.3 Magnetic susceptibility values for an ochre deposit shown with a geological cross

**7.1 Seismic evaluation of surface risks in developed areas overlying former mine** 

Closed or abandoned coal mines are often monitored using seismic surveys, which can locate and characterize workings a few meters below the surface, to depths of up to 60 m. Such workings may collapse or destabilize their surroundings, and thus pose a significant risk to surface activities. The stress field surrounding an abandoned working depends on the initial stress affecting the subsurface structure, which is in turn related to the depth and geometry of the working. The stress field can be evident as both continuous and discontinuous deformation (Mange & Kochonov, 1994). Continuous deformation occurs as elastic and plastic strain in layers, which deform but maintain their continuity. Discontinuous deformation manifests as fractures and displacement of bedrock, evident in

**6.4 Magnetic susceptibility studies of ochra deposit** 

900 x 10-5 SI units (Kotlarczyk & Ratajczak, 2002) (Fig. 6.3).

section of the deposit (Kotlarczyk & Ratajczak, 2002).

**sites** 

**7. Environmental seismic investigations (Jerzy Dec)** 

small scale faulting and other discontinuities (Goodman & Shi, 1985). Discontinuous deformational features pose the greatest risk to the overlying surface. In solid rocks, seismic surveys can identify fractured zones and potential fault planes as velocity changes in seismic data. A seismic refraction survey of shallow workings at a former mine site within Carboniferous rocks, revealed differing behaviour and degrees of deformation (Table 7.1).


Table 7.1 Seismic categories of Carboniferous rocks

## **7.1.1 Seismic surveys in developed areas**

Standard seismic surveys often cannot be conducted at construction sites in developed areas due to their disruptiveness, risks to nearby inhabitants, and accessibility of the area surrounding the site.

In these cases, the uppermost body of bedrock beneath the surface can be imaged by seismic refraction tomography (Fig. 7.1). This technique interprets changes in wave velocity during P-wave propagation, in this example, through a near-surface Quaternary deposit (AB and CD in Fig. 7.1; velocity = V0) and underlying bedrock (BC in Fig. 7.1; velocity = V1). P-wave velocity for the Quaternary deposit was much smaller than that of the Carboniferous rocks (V0 < V1), which resulted in a relatively small critical angle of refraction (i < 10o). The near surface Quaternary deposit was only a few meters thick. AB was thus insignificant, and the distance BC was assumed to equal that of AD. The errors for estimating unit thickness and their associated velocity changes were based on critical angle, i, and refractor depth, h.

Given an AD thickness of 100 m, refractor depth h of 2 - 8m, and i equal to 5 - 15o, errors were relatively small, ranging from 0.3 – 4.3 %. In the model described above, the initial arrival time for the refracted signal *TABCD* (traversing ABCD), could be approximated as the time necessary to traverse BC (Dec, 2004). The study used a seismic source located at point B, and receiver at point C, which reduced the initial arrival time to *TAB*, or time necessary to traverse distance AB. Arrival times were used to construct a tomographic image of the Pwave velocity distribution for rocks beneath the Carboniferous surface. Velocity values were

Geophysics in Near-Surface Investigations 71

Quaternary deposits that overly the Carboniferous unit include several meters of loamy sands and clays. The Carboniferous unit is comprised of sandstone, mudstone, shale and coal seams. The contact between the two units is an angular unconformity (lower units dip ~ 6 - 8o) with a pronounced erosional surface. The survey showed numerous locations where the rock mass was continuous and not fractured, as indicated by velocities of VP > 1900 m/s (Table 7.1). Most of the subsurface in the study area however was strongly fractured and appeared as a discontinuous refraction boundary (700>VP>1800 m/s). A representative location of the study area was the Katowice-D district (Fig. 7.2), where the Carboniferous unit exists in each of the different physical states or categories as described in Table 7.1.

Velocity, VP, from profile 15 was equal to around ~2200 to 1100 m/s (Fig. 7.3), and the state of the Carboniferous surface (B-C2) was categorized as significantly fractured near the edges of the profile. The adjacent profile 16 had VP=900 m/s, and rock structure category D, having no contiguous structure. Data from boreholes drilled in the survey area verified results of the seismic survey. Concrete was then injected into the top parts of the Carboniferous unit through the boreholes. The concrete was readily absorbed by the

The 160 x 70 m investigation area was located between two streets, and was thus divided into two sections (Fig. 7.2). Measurements were performed using a TERRALOC MK6 seismic system. Accelerated weight-drop source EWG-III and L-40B geophones (100 Hz) were used to record seismic sources. Descriptive results for the structure and physical state of the Carboniferous rocks are presented as categories in (Table 7.1). Figure 7.2 illustrates that the subsurface structure was unstable throughout most of study area (category B-D). Two different zones of the study area were distinguished based on the physical state of the subsurface (Fig. 7.2). The left-hand zone was characterized by fractured rock (category B-C2) and the right-hand zone was characterized by rocks that had been pulverized (category D). A contour outlining material with velocity VP = 900 m/s, separates the zone with a fractured foundation. A low velocity anomaly (VP < 700 m/s) is probably related to old workings (cavities) which connect to the ground surface through a system of fractures. Such zones are at risk of infiltration, suffusion and ground collapse. In conclusion, refraction tomography was effective in determining rock properties beneath a developed area, where traditional seismic methods could not be used. Fracture zones overlying old mine workings were evident as strong velocity anomalies, whose magnitude could be correlated to the degree of rock fracturing. Seismic refraction methods are thus suitable for evaluating surface risks of

Fig. 7.3 Cross-sections along refraction profile 15, Katowice-D, Poland

fractures and cavities.

abandoned mines in developed areas.

directly correlated to the composition and physical state of the subsurface material. Fractures and void space for example were evident as areas of reduced velocity.

Fig. 7.1 Seismic refraction tomography: the study model

#### **7.1.2 Evaluating ground surface stability risks**

A survey carried out at the K-K coal mine in Katowice, Poland was used to image the continuity and stability of a subsurface Carboniferous unit, and the risk it posed to surface activities. The survey imaged old shallow (20-50 m) workings dating from the late 19th to the 20th century, which were located in upper parts of the Carboniferous unit. The workings have caused extensive fracturing of the unit (VP <1200 m/s) and cavity migration towards the surface. Fractures from the working significantly increase the risk of ground surface deformation and failure. Due to urban development in the area, a standard seismic profile was not possible; refraction seismic tomography was used instead.

Fig. 7.2 Location and results (velocity distribution) of a tomographic study conducted in a developed area, Katowice-D, Poland; colours indicate velocity values [km/s] from the Carboniferous surface

directly correlated to the composition and physical state of the subsurface material.

A survey carried out at the K-K coal mine in Katowice, Poland was used to image the continuity and stability of a subsurface Carboniferous unit, and the risk it posed to surface activities. The survey imaged old shallow (20-50 m) workings dating from the late 19th to the 20th century, which were located in upper parts of the Carboniferous unit. The workings have caused extensive fracturing of the unit (VP <1200 m/s) and cavity migration towards the surface. Fractures from the working significantly increase the risk of ground surface deformation and failure. Due to urban development in the area, a standard seismic profile

Fig. 7.2 Location and results (velocity distribution) of a tomographic study conducted in a developed area, Katowice-D, Poland; colours indicate velocity values [km/s] from the

Fractures and void space for example were evident as areas of reduced velocity.

Fig. 7.1 Seismic refraction tomography: the study model

was not possible; refraction seismic tomography was used instead.

**7.1.2 Evaluating ground surface stability risks** 

Carboniferous surface

Quaternary deposits that overly the Carboniferous unit include several meters of loamy sands and clays. The Carboniferous unit is comprised of sandstone, mudstone, shale and coal seams. The contact between the two units is an angular unconformity (lower units dip ~ 6 - 8o) with a pronounced erosional surface. The survey showed numerous locations where the rock mass was continuous and not fractured, as indicated by velocities of VP > 1900 m/s (Table 7.1). Most of the subsurface in the study area however was strongly fractured and appeared as a discontinuous refraction boundary (700>VP>1800 m/s). A representative location of the study area was the Katowice-D district (Fig. 7.2), where the Carboniferous unit exists in each of the different physical states or categories as described in Table 7.1.

Velocity, VP, from profile 15 was equal to around ~2200 to 1100 m/s (Fig. 7.3), and the state of the Carboniferous surface (B-C2) was categorized as significantly fractured near the edges of the profile. The adjacent profile 16 had VP=900 m/s, and rock structure category D, having no contiguous structure. Data from boreholes drilled in the survey area verified results of the seismic survey. Concrete was then injected into the top parts of the Carboniferous unit through the boreholes. The concrete was readily absorbed by the fractures and cavities.

Fig. 7.3 Cross-sections along refraction profile 15, Katowice-D, Poland

The 160 x 70 m investigation area was located between two streets, and was thus divided into two sections (Fig. 7.2). Measurements were performed using a TERRALOC MK6 seismic system. Accelerated weight-drop source EWG-III and L-40B geophones (100 Hz) were used to record seismic sources. Descriptive results for the structure and physical state of the Carboniferous rocks are presented as categories in (Table 7.1). Figure 7.2 illustrates that the subsurface structure was unstable throughout most of study area (category B-D). Two different zones of the study area were distinguished based on the physical state of the subsurface (Fig. 7.2). The left-hand zone was characterized by fractured rock (category B-C2) and the right-hand zone was characterized by rocks that had been pulverized (category D). A contour outlining material with velocity VP = 900 m/s, separates the zone with a fractured foundation. A low velocity anomaly (VP < 700 m/s) is probably related to old workings (cavities) which connect to the ground surface through a system of fractures. Such zones are at risk of infiltration, suffusion and ground collapse. In conclusion, refraction tomography was effective in determining rock properties beneath a developed area, where traditional seismic methods could not be used. Fracture zones overlying old mine workings were evident as strong velocity anomalies, whose magnitude could be correlated to the degree of rock fracturing. Seismic refraction methods are thus suitable for evaluating surface risks of abandoned mines in developed areas.

Geophysics in Near-Surface Investigations 73

equipment. This procedure combines signals recorded from the same seismic source and at the same geophone location (Upadhyay, 2004). As many as 20 stacked signals were required to generate high resolution records for some areas. A model of the deposit and the seismic profile obtained from the 24 channel end-off spread are shown in figure 7.4. The following data processing procedures were applied to enhance profile resolution: field static correction, spherical divergence (spreading) correction, refraction static based on the first break picking, surface consistent scaling, surface consistent deconvolution, shaping filter based on the refraction wavelet, Ormsby filter 40/60-150/200 Hz, iterative velocity analysis, NMO, residual static correction, and CMP stacking. Several profiles, referred to as time sections, were recorded throughout the extraction process to monitor temporal changes

Fig. 7.4 Deposit model and seismic raw record; 1-top of the deposit, 4 sand horizon in the

Fig. 7.5 Osiek mine, location of seismic profiles; 1-unexploited zone, 2-producing wells,

3-wells prepared for exploitation, 4-wells with increasing temperature

associated with various phases of extraction.

overburden

#### **7.2 High resolution seismic techniques for investigating the economic sulphur deposits**

Exploitation of sulphur deposits using the underground 'melting' method (also referred to as well mining, or hydrodynamic process mining, or Frash method) poses significant environmental risks. A seismic investigation of the Osiek sulphur mine (Staszow district, Poland) illustrates the usefulness of seismic surveys in managing these risks. Sulphur is extracted at Osiek by pumping hot extraction fluids into the subsurface sulphur deposit. Melted, liquid sulphur is then brought to the surface through a series of production wells. Although the Osiek mine is the only operation in the world to use this technique, extraction at Osiek as well as similar oil and gas operations, put soils, ground water, and air quality at risk. Mine operators can manage risks by closely monitoring the melting process, and incorporating real time geophysical data into risk management decisions. Seismic surveys can specifically be used to monitor the spatial range of the melting zone (Dec, 2010) (i.e. its rate of expansion), and potential subsidence of the overburden (materials overlying the deposit). A seismic investigation of the Osiek mine site illustrates the application described above. Fluid removal of sulphur at the Osiek mine decreases the mechanical strength of porous limestone host rock, which can fail under the pressure of the overburden. Melting zones are not only at risk of failure, but melting may also cause changes in elastic properties of the areas surrounding the extraction zone. Such changes can induce subsidence and other disruptions to the subsurface stress field. High-resolution reflected seismic profiles provide accurate images of the structure of the deposit, its overburden, and signs of deformation in the surroundings. Seismic methods were also used to monitor the melting front and to watch for potential initiation of subsidence troughs, whose risks can be better managed by early detection (Al-Rawahy & Goulty, 1995). Changes in subsurface structure imparted by the Osiek melting method were detected by seismic surveys as early as a month after melting operations began at a particular well site. The extraction fluids used to control the size and expansion rate of the melt zone can also be used to manage subsidence troughs. Under these circumstances, pressure and flow direction of the extraction fluids were used to shape the subsidence trough in a way that prevented the destruction of adjacent wells located along the trough margins.

#### **7.2.1 Field methods for the Osiek mine seismic survey**

Seismic measurements imaged the Miocene unit being exploited by the mining operation and underlying sub-Tertiary strata, at depths of 20 – 200 m. The survey was conducted along lines defined by the positions of overlying production wells (Fig. 7.5). A short spread, high frequency seismic source (Elastic Wave Generator EWG-III, 250 kg accelerated weightdrop), and high frequency geophones (100 Hz) were used to obtain high vertical and horizontal resolution (Brouwer & Helbig, 1998). In order to avoid errors caused by heterogeneities in the subsurface, a split spread with a 100 m interval was applied. Geophones were positioned at 5 m intervals (trace spacing). Near offset was equal to 50 m and far offset was equal to 165 m. The 5 m trace spacing and 1 m shot interval gave CMP resolution equal to 2.5 m and a 24 fold CMP gather.

The majority of the profile was collected using the full 24 fold CMP gather. For parts of the study area that were difficult to access, a 24 channel end-off spread with a 12 fold CMP gather was used. Vertical stacking of signals helped reduce ambient noise caused by mine

Exploitation of sulphur deposits using the underground 'melting' method (also referred to as well mining, or hydrodynamic process mining, or Frash method) poses significant environmental risks. A seismic investigation of the Osiek sulphur mine (Staszow district, Poland) illustrates the usefulness of seismic surveys in managing these risks. Sulphur is extracted at Osiek by pumping hot extraction fluids into the subsurface sulphur deposit. Melted, liquid sulphur is then brought to the surface through a series of production wells. Although the Osiek mine is the only operation in the world to use this technique, extraction at Osiek as well as similar oil and gas operations, put soils, ground water, and air quality at risk. Mine operators can manage risks by closely monitoring the melting process, and incorporating real time geophysical data into risk management decisions. Seismic surveys can specifically be used to monitor the spatial range of the melting zone (Dec, 2010) (i.e. its rate of expansion), and potential subsidence of the overburden (materials overlying the deposit). A seismic investigation of the Osiek mine site illustrates the application described above. Fluid removal of sulphur at the Osiek mine decreases the mechanical strength of porous limestone host rock, which can fail under the pressure of the overburden. Melting zones are not only at risk of failure, but melting may also cause changes in elastic properties of the areas surrounding the extraction zone. Such changes can induce subsidence and other disruptions to the subsurface stress field. High-resolution reflected seismic profiles provide accurate images of the structure of the deposit, its overburden, and signs of deformation in the surroundings. Seismic methods were also used to monitor the melting front and to watch for potential initiation of subsidence troughs, whose risks can be better managed by early detection (Al-Rawahy & Goulty, 1995). Changes in subsurface structure imparted by the Osiek melting method were detected by seismic surveys as early as a month after melting operations began at a particular well site. The extraction fluids used to control the size and expansion rate of the melt zone can also be used to manage subsidence troughs. Under these circumstances, pressure and flow direction of the extraction fluids were used to shape the subsidence trough in a way that prevented the destruction of adjacent wells

Seismic measurements imaged the Miocene unit being exploited by the mining operation and underlying sub-Tertiary strata, at depths of 20 – 200 m. The survey was conducted along lines defined by the positions of overlying production wells (Fig. 7.5). A short spread, high frequency seismic source (Elastic Wave Generator EWG-III, 250 kg accelerated weightdrop), and high frequency geophones (100 Hz) were used to obtain high vertical and horizontal resolution (Brouwer & Helbig, 1998). In order to avoid errors caused by heterogeneities in the subsurface, a split spread with a 100 m interval was applied. Geophones were positioned at 5 m intervals (trace spacing). Near offset was equal to 50 m and far offset was equal to 165 m. The 5 m trace spacing and 1 m shot interval gave CMP

The majority of the profile was collected using the full 24 fold CMP gather. For parts of the study area that were difficult to access, a 24 channel end-off spread with a 12 fold CMP gather was used. Vertical stacking of signals helped reduce ambient noise caused by mine

**7.2 High resolution seismic techniques for investigating the economic sulphur** 

**deposits** 

located along the trough margins.

**7.2.1 Field methods for the Osiek mine seismic survey** 

resolution equal to 2.5 m and a 24 fold CMP gather.

equipment. This procedure combines signals recorded from the same seismic source and at the same geophone location (Upadhyay, 2004). As many as 20 stacked signals were required to generate high resolution records for some areas. A model of the deposit and the seismic profile obtained from the 24 channel end-off spread are shown in figure 7.4. The following data processing procedures were applied to enhance profile resolution: field static correction, spherical divergence (spreading) correction, refraction static based on the first break picking, surface consistent scaling, surface consistent deconvolution, shaping filter based on the refraction wavelet, Ormsby filter 40/60-150/200 Hz, iterative velocity analysis, NMO, residual static correction, and CMP stacking. Several profiles, referred to as time sections, were recorded throughout the extraction process to monitor temporal changes associated with various phases of extraction.

Fig. 7.4 Deposit model and seismic raw record; 1-top of the deposit, 4 sand horizon in the overburden

Fig. 7.5 Osiek mine, location of seismic profiles; 1-unexploited zone, 2-producing wells, 3-wells prepared for exploitation, 4-wells with increasing temperature

Geophysics in Near-Surface Investigations 75

Fig. 7.7 a) Seismic image of overburden layers and associated deformation zone;

labels the same as those in Fig. 7.6

**7.2.3 Studies of water flow in the Osiek subsurface** 

fluid infiltration into the overburden strata.

b) Subsidence effects around exploitation zone and the area of subsidence trough; numerical

Infiltration of fluids into the overburden can change the elastic properties of the materials therein, and thus cause the disappearance of natural reflectors observed in previous surveys. For example, fluid migration into arenaceous horizons of the overburden caused the disappearance of a reflector, demonstrating the physical response of subsurface materials to infiltration. The water infiltration occurred within a zone of the overburden that allowed investigators to identify the source of the infiltration. Successful location of the leak allowed for repair of the failed component. Subsurface monitoring of infiltration is recommended if a given recharge point (e.g. the C-45 recharge well) is located a significant distance from its discharge (eruption) zone (Fig. 7.8). The unexploited part of the deposit is located to the left and the melting-induced destruction zone is located on the right side of figure 7.8. Strong subsidence above the exploited area caused failure of the C-45 well and

The disappearance of seismic reflectors in profiles affected by the infiltration event demonstrates that the zone of infiltration can be precisely identified (Fig. 7.9). Near the F-71 well for example, all seismic reflectors found within the overburden strata disappeared. Fluid infiltration of the arenaceous horizon at 95 m depth was also visible. Post infiltration seismic profiles also indicated that subsidence associated with exploitation of this particular area may facilitate water migration not only within the arenaceous horizon, but also through horizons of the overburden. The broad infiltrated zone around the C-45 well narrowed into an elongated flow path directed to the F-71 well, a site of intensive surface discharge. The system of parallel

seismic sections allowed precise imaging of the direction and width of the flow path.

#### **7.2.2 Influence of melting on the overburden**

Time sections obtained prior to the initiation of melting revealed a number of continuous reflectors in the overburden. Profiles collected during extraction show distinct changes in the top surface of the deposit (above the producing formation), which usually appeared as discontinuities in previously continuous reflectors. These changes were also evident as changes in amplitudes and even disappearance of certain reflectors.

Fig. 7.6 Seismic sections of the sulphur –bearing horizon obtained before (a) and during exploitation (b); 1-top of the sulphur horizon, 2-unexploited zone, 3-exploitation zone, 4 horizons in overburden, 5-deformation zone

Anomalies were also evident in seismic images of the overburden material. These anomalies (50-100 ms) appeared as discontinuities in specific reflectors, indicating subsidence of the formation. Time sections of profile 15 collected prior to initiation of melting, and 3 months after initiation of melting, showed changes in the subsurface (Fig. 7.6). Asymmetric subsidence resulted from destruction and collapse of the overburden strata, which was inturn related to the expansion rate and direction of the melting front. Intuitively, observed deformation was more severe in areas having more closely spaced production wells (Fig. .7a). Deformation of the overburden also appeared to correlate with the boundary between exploited and unexploited zones of the deposit. The most intensive deformation and overburden subsidence (broken continuity of strata; to the right of 1609 well), corresponded to the boundary of an unexploited part of the deposit. The 1609 well was destroyed under similar conditions, by stresses generated in the deformation zone. If deformation-induced stresses generated along planes within the overburden formation attain critical values, production wells within the subsidence zone can be subjected to shearing and ultimate destruction. This scenario (Fig. 7.7b) is evident in the time section profile presented in figure 7.7a.

Time sections obtained prior to the initiation of melting revealed a number of continuous reflectors in the overburden. Profiles collected during extraction show distinct changes in the top surface of the deposit (above the producing formation), which usually appeared as discontinuities in previously continuous reflectors. These changes were also evident as

Fig. 7.6 Seismic sections of the sulphur –bearing horizon obtained before (a) and during exploitation (b); 1-top of the sulphur horizon, 2-unexploited zone, 3-exploitation zone, 4-

Anomalies were also evident in seismic images of the overburden material. These anomalies (50-100 ms) appeared as discontinuities in specific reflectors, indicating subsidence of the formation. Time sections of profile 15 collected prior to initiation of melting, and 3 months after initiation of melting, showed changes in the subsurface (Fig. 7.6). Asymmetric subsidence resulted from destruction and collapse of the overburden strata, which was inturn related to the expansion rate and direction of the melting front. Intuitively, observed deformation was more severe in areas having more closely spaced production wells (Fig. .7a). Deformation of the overburden also appeared to correlate with the boundary between exploited and unexploited zones of the deposit. The most intensive deformation and overburden subsidence (broken continuity of strata; to the right of 1609 well), corresponded to the boundary of an unexploited part of the deposit. The 1609 well was destroyed under similar conditions, by stresses generated in the deformation zone. If deformation-induced stresses generated along planes within the overburden formation attain critical values, production wells within the subsidence zone can be subjected to shearing and ultimate destruction. This scenario (Fig. 7.7b) is evident in the time section

**7.2.2 Influence of melting on the overburden** 

horizons in overburden, 5-deformation zone

profile presented in figure 7.7a.

changes in amplitudes and even disappearance of certain reflectors.

Fig. 7.7 a) Seismic image of overburden layers and associated deformation zone; b) Subsidence effects around exploitation zone and the area of subsidence trough; numerical labels the same as those in Fig. 7.6

#### **7.2.3 Studies of water flow in the Osiek subsurface**

Infiltration of fluids into the overburden can change the elastic properties of the materials therein, and thus cause the disappearance of natural reflectors observed in previous surveys. For example, fluid migration into arenaceous horizons of the overburden caused the disappearance of a reflector, demonstrating the physical response of subsurface materials to infiltration. The water infiltration occurred within a zone of the overburden that allowed investigators to identify the source of the infiltration. Successful location of the leak allowed for repair of the failed component. Subsurface monitoring of infiltration is recommended if a given recharge point (e.g. the C-45 recharge well) is located a significant distance from its discharge (eruption) zone (Fig. 7.8). The unexploited part of the deposit is located to the left and the melting-induced destruction zone is located on the right side of figure 7.8. Strong subsidence above the exploited area caused failure of the C-45 well and fluid infiltration into the overburden strata.

The disappearance of seismic reflectors in profiles affected by the infiltration event demonstrates that the zone of infiltration can be precisely identified (Fig. 7.9). Near the F-71 well for example, all seismic reflectors found within the overburden strata disappeared. Fluid infiltration of the arenaceous horizon at 95 m depth was also visible. Post infiltration seismic profiles also indicated that subsidence associated with exploitation of this particular area may facilitate water migration not only within the arenaceous horizon, but also through horizons of the overburden. The broad infiltrated zone around the C-45 well narrowed into an elongated flow path directed to the F-71 well, a site of intensive surface discharge. The system of parallel seismic sections allowed precise imaging of the direction and width of the flow path.

Geophysics in Near-Surface Investigations 77

Geoelectric profiles were strongly influenced by the water saturated overburden, and were consistent with seismic profiles. The GPR survey showed several areas of the overburden where shallow layers of differing inclination highlighted the landslide's progression. Large blocks of rocks within the slide probably caused the irregular movement observed in the overall landslide. Visible variations in bedrock morphology (Fig. 7.10) provided evidence that the landslide plane formed within flysch strata. Geophysical investigations enabled determination of the depth and horizontal range of the landslide. The integrated interpretation was used to select sites for geotechnical boreholes, and forecast the landslide's future stability. These interpretations suggested that even slight shifts in geotechnical features could reactivate the landslide and jeopardize a nearby road. The risk is a natural consequence of the geomorphology of the study area, where fine-grained units within the

The examples presented above demonstrate how geophysical methods can be used in near surface applications such as resource development, engineering, archaeology, mitigation, remediation and environmental protection. In all cases, the available geophysical techniques were applied using modified methodologies, field measurements and processing strategies. Contemporary catastrophic events and Earth processes in macro scale force application of current geophysical methods in recognition of their results influencing the environment. Heavy rains and floods in Poland in the last ten years gave rise to intensive works in identification and recognizing and mitigation of landslides. Results of geophysical surveys including GPR and shallow seismic refraction were integrated for investigation of landslides. The presented investigations enabled determination of the depth and horizontal range of the landslide. The integrated interpretation showed sites for geotechnical boreholes, and forecasted the landslide's future stability. Geophysical surveys were also useful in preparing the strategy for landslide risk mitigation and helped in the detailed understanding of the internal structure of the landslides, especially

Fig. 7.10 Landslide - seismic refraction depth section.

subsurface can easily become natural slide surfaces.

**8. Conclusions** 

the slide surface.

Fig. 7.8 Seismic section in the zone around the failed well; 1-top of the bed, 2-unexploited zone, 3-exploitation zone, 4-horizons in the overburden, 5-deformation zone, 6-recharge zone.

Fig. 7.9 Seismic profile of a zone affected by extraction fluid migration from the overburden to the surface; labels 1 through 6 are the same as in Fig. 7.8; 7-flow and eruption zone.

#### **7.3 Evaluating landslide stability from seismic surveys**

The current geophysical methods for investigating landslides assume that each landslide must be analyzed separately. Various geophysical methods, including seismic refraction, geoelectrical profiles, and GPR can be integrated in geophysical analysis of landslides. Integrated methods were used to study a landslide that occurred in the Krynica subunit of the Magura unit, Outer Carpathians. The landslide covered approximately 3 hectares, and was located on a slope that descends to a river. The slide surface was initiated in thin, weathered subsurface rocks, and in outcrops of the Eocene Piwniczna sandstones (the sandstone member of the Magura unit). The Piwniczna sandstone is a massive, thickly bedded (1.5 - 6 m) unit with conglomeratic horizons. Fractures are rare in the unit. Within the Piwniczna, packets of medium to fine grained, rhythmically bedded flysch occur. These beds are generally several meters thick. Thin-bedded sandstones and marly shales were also observed in the eastern part of landslide. Field measurements of the landslide were difficult due to the steep gradient of the slope, dense vegetation, and damage to slide surface caused by recent minor reactivation. Nevertheless, eight seismic refraction profiles, three geoelectrical profiles, and six GPR profiles were collected from accessible areas. The seismic boundary separating the overburden (VP ~ 1000 - 1400 m/s) and basement (VP ~ 2200 - 2500 m/s) were distinctly visible on seismic records (Fig. 7.10).

Fig. 7.8 Seismic section in the zone around the failed well; 1-top of the bed, 2-unexploited zone, 3-exploitation zone, 4-horizons in the overburden, 5-deformation zone, 6-recharge

Fig. 7.9 Seismic profile of a zone affected by extraction fluid migration from the overburden to the surface; labels 1 through 6 are the same as in Fig. 7.8; 7-flow and eruption zone.

The current geophysical methods for investigating landslides assume that each landslide must be analyzed separately. Various geophysical methods, including seismic refraction, geoelectrical profiles, and GPR can be integrated in geophysical analysis of landslides. Integrated methods were used to study a landslide that occurred in the Krynica subunit of the Magura unit, Outer Carpathians. The landslide covered approximately 3 hectares, and was located on a slope that descends to a river. The slide surface was initiated in thin, weathered subsurface rocks, and in outcrops of the Eocene Piwniczna sandstones (the sandstone member of the Magura unit). The Piwniczna sandstone is a massive, thickly bedded (1.5 - 6 m) unit with conglomeratic horizons. Fractures are rare in the unit. Within the Piwniczna, packets of medium to fine grained, rhythmically bedded flysch occur. These beds are generally several meters thick. Thin-bedded sandstones and marly shales were also observed in the eastern part of landslide. Field measurements of the landslide were difficult due to the steep gradient of the slope, dense vegetation, and damage to slide surface caused by recent minor reactivation. Nevertheless, eight seismic refraction profiles, three geoelectrical profiles, and six GPR profiles were collected from accessible areas. The seismic boundary separating the overburden (VP ~ 1000 - 1400 m/s) and basement (VP ~ 2200 - 2500

**7.3 Evaluating landslide stability from seismic surveys** 

m/s) were distinctly visible on seismic records (Fig. 7.10).

zone.

Fig. 7.10 Landslide - seismic refraction depth section.

Geoelectric profiles were strongly influenced by the water saturated overburden, and were consistent with seismic profiles. The GPR survey showed several areas of the overburden where shallow layers of differing inclination highlighted the landslide's progression. Large blocks of rocks within the slide probably caused the irregular movement observed in the overall landslide. Visible variations in bedrock morphology (Fig. 7.10) provided evidence that the landslide plane formed within flysch strata. Geophysical investigations enabled determination of the depth and horizontal range of the landslide. The integrated interpretation was used to select sites for geotechnical boreholes, and forecast the landslide's future stability. These interpretations suggested that even slight shifts in geotechnical features could reactivate the landslide and jeopardize a nearby road. The risk is a natural consequence of the geomorphology of the study area, where fine-grained units within the subsurface can easily become natural slide surfaces.

## **8. Conclusions**

The examples presented above demonstrate how geophysical methods can be used in near surface applications such as resource development, engineering, archaeology, mitigation, remediation and environmental protection. In all cases, the available geophysical techniques were applied using modified methodologies, field measurements and processing strategies.

Contemporary catastrophic events and Earth processes in macro scale force application of current geophysical methods in recognition of their results influencing the environment. Heavy rains and floods in Poland in the last ten years gave rise to intensive works in identification and recognizing and mitigation of landslides. Results of geophysical surveys including GPR and shallow seismic refraction were integrated for investigation of landslides. The presented investigations enabled determination of the depth and horizontal range of the landslide. The integrated interpretation showed sites for geotechnical boreholes, and forecasted the landslide's future stability. Geophysical surveys were also useful in preparing the strategy for landslide risk mitigation and helped in the detailed understanding of the internal structure of the landslides, especially the slide surface.

Geophysics in Near-Surface Investigations 79

Dec, J. (2010). High Resolution Seismic Investigations for the Determination of Water Flow

Gołębiowski, T., Tomecka-Suchoń, S., Marcak, H. & Żogała, B. (2010a). Aiding of the GPR

Gołębiowski, T., Marcak, H., Tomecka-Suchoń, S., Zdechlik, R., Zuberek, W. (2010b). Use of

Goodman, E., E. & SHI G., H. (1985). *Block Theory and its Application to Rock Engineering*.

Halliburton (1991). *Log Interpretation Charts,* Halliburton Logging Services, Inc., Houston,

Hearst, J., R., Nelson, P., H. & Paillet, F., L. (2000). *Well Logging for Physical Properties.* 

Jarzyna, J., Bała, M. & Cichy, A. (2010). Elastic parameters of rocks from well logging in

Jarzyna, J., Puskarczyk, E., Bała, M. & Papiernik, B. (2009). Variability of the Rotliegend

Kobranova, V., N. (1986). (Eng. Translation 1989). *Petrophysics*. MIR Publish. Moscow, Springer-Verlag, ISBN 3-540-51524-0, Berlin, Heidelberg, New York. Kotlarczyk, J. & Ratajczak, T. (2002). *Carpathian ochre from Czerwonki Hermanowskie near* 

Liu Y.X., (2007). Evaluation and extraction of weak gravity and magnetic anomalies*. Applied* 

Madej J., Jakiel K. & Porzucek S. (2001). Microgravimetric assesment of possible surface

Mange, C. & Kochonov M. (1997). Anisotropic Material with Interacting Arbitrarily

Marcak, H. & Tomecka-Suchoń S. (2010). Properties of Georadar Signals Used for an

Nguyen, van, G., Ziętek, J., Nguyen, Ba, D., Karczewski, J. & Gołębiowski, T. (2005). Study

Near Surface Sediments. *Acta Geophysica*, Vol. 58, No 1, pp. 34-48.

*A Handbook for Geophysicists, Geologists and Engineers*. Wiley & Sons, ISBN 0-471-

Sandstones in the Polish Part of the Southern Permian Basin - Permeability and Porosity Relationships. *Annales Societatis Geologorum Poloniae*, Vol. 79, pp. 13-26. Jol , H., M. (2009). *Ground Penetrating Radar Theory and Applications*, Elsevier, ISBN 978-0-444-

*Tyczyn*, Mineral and Energy Economy Research Institute, ISBN 83-89174-60-X,

deformations in some post-mining areas. Mineral Deposits at the Beginning of the 21st Century, Proc. of the Joint Sixth Biennial SGA-SEG Meeting, Krakow, 26-29

Oriented Cracks. Stress Intensity Factor and Crack-Microcrack Interactions.

Estimation of the Mineralization of the Soil Waters, *Archives of Mining Sciences*, Vol.

of Geological Sedimentary Structures of the Mekong River Banks by Ground Penetrating Radar: Forecasting Avulsion-Prone Zones, Acta Geophysica Polonica,

21-25 June 2010.

Texas, USA.

Eds: Kania J., Kmiecik, E., Zuber A.

96305-4, Chichester, England.

53348-7, Amsterdam.

Krakow (in Polish).

55, No. 3, pp. 469-487.

*Geophysics*, Vol. 4, No.4, pp. 288-293.

August 2001, A.A. Balkema Publishers, pp.1035-1038.

*International Journal of Fracture*, Vol. 65, pp. 115-141.

Vol. 53, No. 2, pp. 167-181, ISSN 1895-6572.

Englewood Clifs, N.J., Prentice-Hall, Inc.

Directions During Sulphur Deposits Exploitation. *Acta Geophysica*, Vol. 58, pp. 5-14.

Method by the Other Measurement Techniques for the Liquid Contamination Detection, XIII International Conference on Ground Penetrating Radar, Lecce/

Geophysical Methods for the Assessment of Migration of Contaminants from the Coal-Mining Waste Dumps, Extended Abstracts, abstract id: 317,,XXXVIII IAH Congress, Groundwater Quality Sustainability, Krakow, 12-17 September 2010,

Results of GPR and microgravity surveys turned out also to be highly effective in locating zones of weakness in river embankments. The results of field measurements enabled qualitative determination of the embankment's inner structure and basal materials. It was shown that GPR and microgravity methods were useful in the investigation of embankments and earthen dams that are at risk of water leakage and possible failure.

Well known applications of GPR in archaeology were also presented together with new area oriented to identification of low conductivity and high conductivity contaminations of subsurface formations.

Seismics in tomography mode and refraction mode and microgravity surveys were shown as methods useful in identification and recognition of weak zones in the industrial areas of former and contemporary mining activity. Fracture zones overlying old mine workings were evident as strong velocity anomalies and low density anomalies, whose magnitude could be correlated to the degree of rock fracturing suitable for evaluating surface risks of abandoned mines in developed areas.

Magnetic susceptibility measurements *in situ* and in laboratory were presented as an useful tool in effective monitoring of anthropogenic influences on the environment in the past (historical layers) and in tracking the contemporary pollution by iron compounds. Magnetic anomalies are also the evidence of weathering of ferrimagnetic and antiferromagnetic minerals, or other pedogenic processes.

Universal petrophysical formulas were demonstrated to combine empirical laboratory data with geophysical and well logging measurements and theoretical expectations and parameters forming relationships useful in a consistent, quantitative geophysical interpretation. Exemplary values of parameters for the lithology types were presented in tables and figures to illustrate the level of parameters values frequently found in the near surface formations.

## **9. Acknowledgments**

Authors would like to thank Ms Teresa Staszowska for preparing the figures and text edition.

## **10. References**


Results of GPR and microgravity surveys turned out also to be highly effective in locating zones of weakness in river embankments. The results of field measurements enabled qualitative determination of the embankment's inner structure and basal materials. It was shown that GPR and microgravity methods were useful in the investigation of embankments and earthen dams that are at risk of water leakage and possible failure.

Well known applications of GPR in archaeology were also presented together with new area oriented to identification of low conductivity and high conductivity contaminations of

Seismics in tomography mode and refraction mode and microgravity surveys were shown as methods useful in identification and recognition of weak zones in the industrial areas of former and contemporary mining activity. Fracture zones overlying old mine workings were evident as strong velocity anomalies and low density anomalies, whose magnitude could be correlated to the degree of rock fracturing suitable for evaluating surface risks of

Magnetic susceptibility measurements *in situ* and in laboratory were presented as an useful tool in effective monitoring of anthropogenic influences on the environment in the past (historical layers) and in tracking the contemporary pollution by iron compounds. Magnetic anomalies are also the evidence of weathering of ferrimagnetic and antiferromagnetic

Universal petrophysical formulas were demonstrated to combine empirical laboratory data with geophysical and well logging measurements and theoretical expectations and parameters forming relationships useful in a consistent, quantitative geophysical interpretation. Exemplary values of parameters for the lithology types were presented in tables and figures to illustrate the level of parameters values frequently found in the near

Authors would like to thank Ms Teresa Staszowska for preparing the figures and text

Al-Rawahy, S., Y., S. & Goulty, N., R. (1995). Effect of Mining Subsidence on Seismic

Annan, A., P. (2001). Ground Penetrating Radar. Workshop Notes, Sensors & Software,

Brouwer, J. & Helbig, K. (1998). *Shallow High-Resolution Reflection Seismics*. Elsevier Science

Daniels, D., J. (2004). *Ground Penetrating Radar* - 2nd Ed., The Institution of Electrical

Dec, J. (2004). Seismic Survey to Evaluate the Danger of Ground Surface Damage in Built-Up Terrain in Mining Areas. *Polish Journal of Environmental Studies*, Vol. 13, pp. 70-73.

Velocity Monitored by a Repeated Reflection Profile, *Geophysical Prospecting,* Vol.

subsurface formations.

surface formations.

**10. References** 

edition.

**9. Acknowledgments** 

43, pp. 191-201.

Ltd., Amsterdam.

Engineers, ISBN 978-0-86341-360-5, GB.

Canada.

abandoned mines in developed areas.

minerals, or other pedogenic processes.


**0**

**4**

F.E.M.(Ted) Lilley

*Australia*

**Magnetotelluric Tensor Decomposition: Insights**

The magnetotelluric (MT) method of geophysics exploits the phenomenon of natural electromagnetic induction which takes place at and near the surface of Earth. The purpose is to determine information about the electrical conductivity structure of Earth, upon which the process of electromagnetic induction depends. The MT method has been well described recently in books such as Simpson & Bahr (2005), Gubbins & Herrero-Bervera (2007), and Berdichevsky & Dmitriev (2008). The reader is referred to these for general information about the method and its results. Notable modern extensions of the method are observation at an array of sites simultaneously, and observation on the seafloor (both shallow and deep oceans). In the most simple form of the method, data are observed at a single field site. Typically, three components (north, east and vertically downwards) of the fluctuating magnetic field are observed, and two components (north and east) of the fluctuating electric field. The magnetic field is measured using a variety of instruments such as fluxgates and induction coils. The electric field is measured more simply, between grounded electrodes typically

The natural signals observed cover a frequency band from 0.001 to 1000 Hz. They have a variety of causes, the relative importance of which varies with position on the Earth, especially latitude. Recorded data are transformed to the frequency domain, and interpretation proceeds

The reduction of observed time-series to the frequency domain is thus fundamental to the MT method. In the frequency domain various transfer functions are determined, encapsulating

Another fundamental part of the data-reduction process, the focus of the present chapter, is "rotation" of observed data. By rotation is meant an examination of the observed transfer functions to see how they would vary were the observing axes rotated at the observing site.

The MT tensor hasa2x2 form and is well-suited to analysis by the methods of linear algebra. This suitability was evident early in the development of MT, and it is common for a paper on MT to start with such an analysis. The papers of Eggers (1982), LaTorraca et al. (1986), and Yee & Paulson (1987), for example, specifically proceed with eigenvalue analysis and singular

the response of the observing site to the source fields causing the induction.

Rotation may reveal geologic dimensionality, and strike direction.

value decomposition (SVD), and see also Weaver (1994).

**1. Introduction**

several hundred metres apart.

based on frequency-dependence.

**from Linear Algebra and Mohr Diagrams**

*Research School of Earth Sciences, Australian National University, Canberra*


## **Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams**

F.E.M.(Ted) Lilley *Research School of Earth Sciences, Australian National University, Canberra Australia*

#### **1. Introduction**

80 New Achievements in Geoscience

Rosowiecka, O. & Nawrocki, J. (2010). Assesment of Soils Pollution Extent in Surroundings

Schön, J., K. (2004). *Physical Properties of Rocks. Fundamentals and Principles of Petrophysics*.

Tiab, D., & Donaldson, E.,C. (2004). *Petrophysics. Theory and Practice of Measuring Reservoir* 

Wojas, A. (2009). Magnetic Susceptibility of Soils, Including Iron Oxides of Anthropogenic

Wójcicki, A. (1993). Approximation of the Gravity Attraction Caused by the Terrain Relief

Francisco, Singapore, Sidney, Tokyo, Gulf Professional Publishing.

Scientists, ISBN 978-83-926896-1-4, Poznań, Poland, May 2009.

No.1, (December 2009), pp.185-194, ISSN 0039-3169.

Upadhyay, S., K. (2004). *Seismic Reflection Processing*, Springer, Berlin.

24.

S., (Eds), v. 18, Elsevier, ISBN 0-08-044346-X, Oxford, UK.

of Ironworks Based on Magnetic Analysis. *Studia Geophysica et Geodaetica*, Vol. 54,

Handbook of Geophysical Exploration, Seismic Exploration, Helbig K. and Treitel

*Rock and Fluid Transport Properties*. (sec. ed.), Elsevier, ISBN 0-7506-7711-2, Amsterdam, Boston, Heidelberg, London, New York, Oxford, Paris, San Diego, San

and Natural Origin - Measurements Using the Bartington Instrument, Proceedings of InterTech - II International Interdisciplinary Technical Conference of Young

Forms Using a Polyhedron Method. *Acta Geophysica Polonica*, Vol. 41, No. 3, pp.1-

The magnetotelluric (MT) method of geophysics exploits the phenomenon of natural electromagnetic induction which takes place at and near the surface of Earth. The purpose is to determine information about the electrical conductivity structure of Earth, upon which the process of electromagnetic induction depends. The MT method has been well described recently in books such as Simpson & Bahr (2005), Gubbins & Herrero-Bervera (2007), and Berdichevsky & Dmitriev (2008). The reader is referred to these for general information about the method and its results. Notable modern extensions of the method are observation at an array of sites simultaneously, and observation on the seafloor (both shallow and deep oceans).

In the most simple form of the method, data are observed at a single field site. Typically, three components (north, east and vertically downwards) of the fluctuating magnetic field are observed, and two components (north and east) of the fluctuating electric field. The magnetic field is measured using a variety of instruments such as fluxgates and induction coils. The electric field is measured more simply, between grounded electrodes typically several hundred metres apart.

The natural signals observed cover a frequency band from 0.001 to 1000 Hz. They have a variety of causes, the relative importance of which varies with position on the Earth, especially latitude. Recorded data are transformed to the frequency domain, and interpretation proceeds based on frequency-dependence.

The reduction of observed time-series to the frequency domain is thus fundamental to the MT method. In the frequency domain various transfer functions are determined, encapsulating the response of the observing site to the source fields causing the induction.

Another fundamental part of the data-reduction process, the focus of the present chapter, is "rotation" of observed data. By rotation is meant an examination of the observed transfer functions to see how they would vary were the observing axes rotated at the observing site. Rotation may reveal geologic dimensionality, and strike direction.

The MT tensor hasa2x2 form and is well-suited to analysis by the methods of linear algebra. This suitability was evident early in the development of MT, and it is common for a paper on MT to start with such an analysis. The papers of Eggers (1982), LaTorraca et al. (1986), and Yee & Paulson (1987), for example, specifically proceed with eigenvalue analysis and singular value decomposition (SVD), and see also Weaver (1994).

O

Fig. 1. The rotation of MT observing axes clockwise by angle *θ*�

 *Z*� *xx Z*� *xy*

*Z*�

*Z*�

*Z*�

*Z*�

*Z*� *yx Z*� *yy*

with **<sup>R</sup>T**(*θ*) the transpose of **<sup>R</sup>**(*θ*). **<sup>R</sup>T**(*θ*) will sometimes be written **<sup>R</sup>**(−*θ*).

Upon rotation of the horizontal measuring axes clockwise by angle *θ*' as shown in Fig. 1,

*Zxx Zxy Zyx Zyy*

*xx*, *Z*� *xy*; *Z*� *yx*, *Z*�

= **R**(*θ*� ) 

Thus the elements of the second matrix are related to the first by the following equations:

*<sup>C</sup>* = [(*Zxx* <sup>−</sup> *Zyy*)<sup>2</sup> + (*Zxy* <sup>+</sup> *Zyx*)2]

east) to OX' and OY'.

where

and *β* is defined by

**2.1 Rotation of the horizontal axes**

matrix [*Zxx*, *Zxy*; *Zyx*, *Zyy*] changes to [*Z*�

X

θ'

X'

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 83

Y

, from OX and OY (north and

) (5)

<sup>2</sup> /2 (10)

Y'

*yy*] according to

**R**(−*θ*�

*xx* = (*Zxx* + *Zyy*)/2 + *C* sin(2*θ*� + *β*) (6)

*xy* = (*Zxy* − *Zyx*)/2 + *C* cos(2*θ*� + *β*) (7)

*yx* = −(*Zxy* − *Zyx*)/2 + *C* cos(2*θ*� + *β*) (8)

*yy* = (*Zxx* + *Zyy*)/2 − *C* sin(2*θ*� + *β*) (9)

1

tan *β* = (*Zxx* − *Zyy*)/(*Zxy* + *Zyx*) (11)

There appear however to be a number of aspects of this approach as yet unexplored, especially when the in-phase and quadrature parts of the tensor are analysed separately. This chapter now follows such a line of enquiry, investigating the significance of symmetric, antisymmetric and non-symmetric parts, eigenvalue analysis and SVD.

The Mohr diagram representation is found to be a useful way to display some of the results. Cases of 1D, 2D and 3D electrical conductivity structure are examined, including the 3D case where the MT phase is greater than 90◦ or "out of quadrant".

Taking the Smith (1995) treatment of telluric distortion, a particular case examined is that where SVD gives an angle which is frequency independent, implying that it contains geological information. Such angles arise when MT tensor data are nearly singular.

The analysis of in-phase and quadrature parts separately is seen to give further insight into the question of whether, as suspected by Lilley (1998), the determinants of the parts taken separately should always be positive.

#### **2. Notation**

The common representation of a magnetotelluric tensor **Z** is taken,

$$\mathbf{E} = \mathbf{Z}\mathbf{H} \tag{1}$$

of components

$$
\begin{bmatrix} E\_{\chi} \\ E\_{y} \end{bmatrix} = \begin{bmatrix} Z\_{\text{xx}} \ Z\_{\text{xy}} \\ Z\_{\text{yx}} \ Z\_{\text{yy}} \end{bmatrix} \begin{bmatrix} H\_{\chi} \\ H\_{y} \end{bmatrix} \tag{2}
$$

linking observed electric **E** and magnetic **H** fluctuations at an observing site on the surface of Earth.

All quantities are complex functions of frequency *ω*, and in Equations 1 and 2 a time dependence of exp(*iωt*) is understood.

In this chapter the subscripts *p*, *q* will be used to denote in-phase and quadrature parts. For example the complex quantity *Z<sup>σ</sup>* is expressed

$$\mathbf{Z}^{\sigma} = \mathbf{Z}\_p^{\sigma} + \mathbf{i}\mathbf{Z}\_q^{\sigma} \tag{3}$$

Note that adopting a time-dependence of exp(−*iωt*) as recommended by authors such as Stratton (1941) and Hobbs (1992) would change the sign of *Z<sup>σ</sup> <sup>q</sup>* . Such a change may be misinterpreted, especially when distortion has caused the phase to be out of its expected first quadrant.

The subscripts *p*, *q* will also be used to denote quantities which are derived from the in-phase and quadrature parts, respectively, of a complex quantity, but which are themselves not recombined to give a further complex quantity. In Sections 2.1, 3, 5 and 11 below, where derivations apply equally to in-phase and quadrature cases, subscripts *p*, *q* are omitted for simplicity. Such derivations can be taken as applying to in-phase components, with similar derivations possible for quadrature components.

Also in this chapter, for compactness of text, a 2 x 2 matrix such as that for **Z** in Equation 2 will in places be written [*Zxx*, *Zxy*; *Zyx*, *Zyy*]. A rotation matrix **R**(*θ*) will be introduced

$$\mathbf{R}(\theta) = [\cos \theta, \sin \theta; -\sin \theta, \cos \theta] \tag{4}$$

Fig. 1. The rotation of MT observing axes clockwise by angle *θ*� , from OX and OY (north and east) to OX' and OY'.

with **<sup>R</sup>T**(*θ*) the transpose of **<sup>R</sup>**(*θ*). **<sup>R</sup>T**(*θ*) will sometimes be written **<sup>R</sup>**(−*θ*).

#### **2.1 Rotation of the horizontal axes**

Upon rotation of the horizontal measuring axes clockwise by angle *θ*' as shown in Fig. 1, matrix [*Zxx*, *Zxy*; *Zyx*, *Zyy*] changes to [*Z*� *xx*, *Z*� *xy*; *Z*� *yx*, *Z*� *yy*] according to

$$
\begin{bmatrix} Z'\_{xx} Z'\_{xy} \\ Z'\_{yx} Z'\_{yy} \end{bmatrix} = \mathbf{R}(\theta') \begin{bmatrix} Z\_{xx} \ Z\_{xy} \\ Z\_{yx} Z\_{yy} \end{bmatrix} \mathbf{R}(-\theta') \tag{5}
$$

Thus the elements of the second matrix are related to the first by the following equations:

$$Z'\_{\rm xx} = (Z\_{\rm xx} + Z\_{yy})/2 + \mathcal{C}\sin(2\theta' + \beta) \tag{6}$$

$$Z\_{xy}^{l} = (Z\_{xy} - Z\_{yx})/2 + \mathcal{C}\cos(2\theta^{l} + \beta) \tag{7}$$

$$Z\_{yx}^{\prime} = -(Z\_{xy} - Z\_{yx})/2 + \mathbb{C}\cos(2\theta^{\prime} + \beta) \tag{8}$$

$$Z\_{yy}' = (Z\_{xx} + Z\_{yy})/2 - \mathbb{C}\sin(2\theta' + \beta) \tag{9}$$

where

2 Will-be-set-by-IN-TECH

There appear however to be a number of aspects of this approach as yet unexplored, especially when the in-phase and quadrature parts of the tensor are analysed separately. This chapter now follows such a line of enquiry, investigating the significance of symmetric, antisymmetric

The Mohr diagram representation is found to be a useful way to display some of the results. Cases of 1D, 2D and 3D electrical conductivity structure are examined, including the 3D case

Taking the Smith (1995) treatment of telluric distortion, a particular case examined is that where SVD gives an angle which is frequency independent, implying that it contains

The analysis of in-phase and quadrature parts separately is seen to give further insight into the question of whether, as suspected by Lilley (1998), the determinants of the parts taken

> *Zxx Zxy Zyx Zyy*

linking observed electric **E** and magnetic **H** fluctuations at an observing site on the surface of

All quantities are complex functions of frequency *ω*, and in Equations 1 and 2 a time

In this chapter the subscripts *p*, *q* will be used to denote in-phase and quadrature parts. For

Note that adopting a time-dependence of exp(−*iωt*) as recommended by authors such as

misinterpreted, especially when distortion has caused the phase to be out of its expected first

The subscripts *p*, *q* will also be used to denote quantities which are derived from the in-phase and quadrature parts, respectively, of a complex quantity, but which are themselves not recombined to give a further complex quantity. In Sections 2.1, 3, 5 and 11 below, where derivations apply equally to in-phase and quadrature cases, subscripts *p*, *q* are omitted for simplicity. Such derivations can be taken as applying to in-phase components, with similar

Also in this chapter, for compactness of text, a 2 x 2 matrix such as that for **Z** in Equation 2 will in places be written [*Zxx*, *Zxy*; *Zyx*, *Zyy*]. A rotation matrix **R**(*θ*) will be introduced

**R**(*θ*)=[cos *θ*, sin *θ*; − sin *θ*, cos *θ*] (4)

*<sup>p</sup>* + *iZ<sup>σ</sup>*

*Z<sup>σ</sup>* = *Z<sup>σ</sup>*

 *Hx Hy* 

**E** = **ZH** (1)

*<sup>q</sup>* (3)

*<sup>q</sup>* . Such a change may be

(2)

geological information. Such angles arise when MT tensor data are nearly singular.

and non-symmetric parts, eigenvalue analysis and SVD.

separately should always be positive.

dependence of exp(*iωt*) is understood.

example the complex quantity *Z<sup>σ</sup>* is expressed

derivations possible for quadrature components.

**2. Notation**

of components

Earth.

quadrant.

where the MT phase is greater than 90◦ or "out of quadrant".

The common representation of a magnetotelluric tensor **Z** is taken,

 *Ex Ey* = 

Stratton (1941) and Hobbs (1992) would change the sign of *Z<sup>σ</sup>*

$$\mathcal{L} = [(Z\_{\rm xx} - Z\_{yy})^2 + (Z\_{\rm xy} + Z\_{yx})^2]^{\frac{1}{2}}/2\tag{10}$$

and *β* is defined by

$$\tan \beta = (Z\_{xx} - Z\_{yy}) / (Z\_{xy} + Z\_{yx}) \tag{11}$$

Z <sup>L</sup>

Z'xy

μ

(Zxx+Zyy)/2

Z'xx

0

a

Z'xxp

b

Z'xxp

c

observed as [0, *Z*; −*Z*, 0].

**3.3 1D structure**

Zp

<sup>χ</sup>Zp

0 0

0 0

reduces to a pair of points on the horizontal axes, as shown in Fig. 2c.

**4. Invariants of rotation of the measuring axes**

marks observed point (Zxy, Zxx)

C

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 85

β' 2θ'

(Zxy-Zyx)/2

Z'xxq

In-phase Quadrature

Z'xxq

Zp Zq

marks general point (Z'xy, Z'xx) after axes rotation of θ'

<sup>σ</sup> Zq

β β

Z'xyq

σ

Z'yy

0

β

Z'yx

<sup>χ</sup>Zq

Z'xyq

Z'xyp

Z'xyp

Fig. 2. Mohr diagrams for a: The general MT tensor; both in-phase and quadrature parts take this general form. b: The 2D case; radial arms are now parallel, with circle centres on the horizontal axes. The 2D values *Z<sup>σ</sup>* and *Z<sup>χ</sup>* could be interchanged. c: The 1D case, for **Z**

If a 2D case simplifies further to become a 1D case, *Z<sup>σ</sup>* = *Z<sup>χ</sup>* = *Z* say, and the tensor for all rotations has the form [0, *Z*; −*Z*, 0]. The length *C* of the radial arm vanishes, and the diagram

It can be seen from Fig. 2a that the MT observations can be expressed in terms of seven invariants of rotation. From different sets which are possible, and evident from formal analysis (Szarka & Menvielle, 1997; Weaver et al., 2000), this chapter adopts the invariants shown in Fig. 3. Thus the eight values of the MT tensor, all of which generally change upon axes rotation, become seven invariants plus one angle. That angle, *θ*� in Fig. 2a, is the angle which

defines the direction of the measuring axes (for example, with respect to north).

It is also useful to define an auxiliary angle *β*� by

$$\tan \beta' = (Z'\_{\rm xx} - Z'\_{yy}) / (Z'\_{\rm xy} + Z'\_{yx}) \tag{12}$$

with

$$\boldsymbol{\theta}' = (\boldsymbol{\beta}' - \boldsymbol{\beta})/2 \tag{13}$$

and angle *μ* as

$$\tan \mu = (Z\_{yy} + Z\_{xx}) / (Z\_{xy} - Z\_{yx}) \tag{14}$$

also *Z<sup>L</sup>* as

$$\mathbf{Z}^{\rm L} = [(\mathbf{Z}\_{\rm xx} + \mathbf{Z}\_{\rm yy})^2 + (\mathbf{Z}\_{\rm xy} - \mathbf{Z}\_{\rm yx})^2]^{\frac{1}{2}}/2 \tag{15}$$

Then plotting *Z*� *xx* against *Z*� *xy* as the axes are rotated (i.e. *θ*� varies) defines a circle, known (with its axes) as a Mohr diagram. On such a figure axes for *Z*� *yx* and *Z*� *yy* may be included, to display the variation of these components also.

#### **3. The depiction of MT tensors using Mohr diagrams**

The Mohr diagram representation, straightforward for a2x2 matrix, is an informative figure for the student of linear algebra. For the MT case it may be a useful way to display results, as will now be shown for the general case of 3D conductivity structure, and the particular cases of 2D and 1D conductivity structure to which the 3D case simplifies.

In this chapter in-phase and quadrature data will be presented in adjacent figures. There is an appeal in putting both in-phase and quadrature data on the same axes, as demonstrated by Szarka & Menvielle (1997) and Weaver et al. (2000), but when showing a full frequency range separate sets of axes are practical. Also at times it is helpful to add axes for *Z*� *yx* and *Z*� *yy* as on Fig. 2a, and such an addition is not possible for *Z*� *xx* and *Z*� *xy* axes common to both in-phase and quadrature parts.

#### **3.1 The general case; 3D structure**

Mohr diagrams for the general case of the MT tensor are shown in Fig. 2a. Different points on the diagram can be checked to confirm that Equations 6 to 15 for the rotation of axes are obeyed. The diagram is "Type 1" of Lilley (1998, p.1889), and shows how axes for *Z*� *yx* and *Z*� *yy* can be included. The diagram in Fig. 2a is drawn for relatively mild 3D characteristics. Diagrams for strong 3D characteristics are discussed below and shown in an example; however, what limits there may be to extreme 3D behaviour have not yet been fully explored. Some examples of 3D behaviour which appear to be prohibited will be discussed in this chapter.

#### **3.2 2D structure**

When the geologic structure is 2D and the axes are rotated to be along and across geologic strike, the MT tensor can be rotated to have the form [0, *<sup>Z</sup>σ*; <sup>−</sup>*Zχ*, 0], where it is expected both *Z<sup>σ</sup>* and *Z<sup>χ</sup>* are positive (a point discussed in Section 13.1 below). As shown in Fig. 2b, the diagram becomes a pair of circles with origins on the *Z*� *xy* axes, and angles *μ<sup>p</sup>* and *μ<sup>q</sup>* are zero. The in-phase and quadrature radial arms are parallel. *Z<sup>σ</sup>* and *Z<sup>χ</sup>* are known as the *TE* and *TM* modes (or vice-versa).

Fig. 2. Mohr diagrams for a: The general MT tensor; both in-phase and quadrature parts take this general form. b: The 2D case; radial arms are now parallel, with circle centres on the horizontal axes. The 2D values *Z<sup>σ</sup>* and *Z<sup>χ</sup>* could be interchanged. c: The 1D case, for **Z** observed as [0, *Z*; −*Z*, 0].

#### **3.3 1D structure**

4 Will-be-set-by-IN-TECH

*xx* − *Z*�

*<sup>Z</sup><sup>L</sup>* = [(*Zxx* <sup>+</sup> *Zyy*)<sup>2</sup> + (*Zxy* <sup>−</sup> *Zyx*)2]

The Mohr diagram representation, straightforward for a2x2 matrix, is an informative figure for the student of linear algebra. For the MT case it may be a useful way to display results, as will now be shown for the general case of 3D conductivity structure, and the particular cases

In this chapter in-phase and quadrature data will be presented in adjacent figures. There is an appeal in putting both in-phase and quadrature data on the same axes, as demonstrated by Szarka & Menvielle (1997) and Weaver et al. (2000), but when showing a full frequency range

Mohr diagrams for the general case of the MT tensor are shown in Fig. 2a. Different points on the diagram can be checked to confirm that Equations 6 to 15 for the rotation of axes are obeyed. The diagram is "Type 1" of Lilley (1998, p.1889), and shows how axes for

characteristics. Diagrams for strong 3D characteristics are discussed below and shown in an example; however, what limits there may be to extreme 3D behaviour have not yet been fully explored. Some examples of 3D behaviour which appear to be prohibited will be discussed in

When the geologic structure is 2D and the axes are rotated to be along and across geologic strike, the MT tensor can be rotated to have the form [0, *<sup>Z</sup>σ*; <sup>−</sup>*Zχ*, 0], where it is expected both *Z<sup>σ</sup>* and *Z<sup>χ</sup>* are positive (a point discussed in Section 13.1 below). As shown in Fig. 2b, the

The in-phase and quadrature radial arms are parallel. *Z<sup>σ</sup>* and *Z<sup>χ</sup>* are known as the *TE* and

*yy* can be included. The diagram in Fig. 2a is drawn for relatively mild 3D

*xx* and *Z*�

*yy*)/(*Z*�

*xy* + *Z*�

*θ*� = (*β*� − *β*)/2 (13)

*yx* and *Z*�

tan *μ* = (*Zyy* + *Zxx*)/(*Zxy* − *Zyx*) (14)

*xy* as the axes are rotated (i.e. *θ*� varies) defines a circle, known

1

*yx*) (12)

<sup>2</sup> /2 (15)

*yy* may be included, to

*yx* and *Z*�

*xy* axes common to both in-phase

*xy* axes, and angles *μ<sup>p</sup>* and *μ<sup>q</sup>* are zero.

*yy* as on

It is also useful to define an auxiliary angle *β*� by

*xx* against *Z*�

display the variation of these components also.

Fig. 2a, and such an addition is not possible for *Z*�

diagram becomes a pair of circles with origins on the *Z*�

(with its axes) as a Mohr diagram. On such a figure axes for *Z*�

**3. The depiction of MT tensors using Mohr diagrams**

of 2D and 1D conductivity structure to which the 3D case simplifies.

separate sets of axes are practical. Also at times it is helpful to add axes for *Z*�

with

and angle *μ* as

Then plotting *Z*�

and quadrature parts.

*Z*�

*yx* and *Z*�

this chapter.

**3.2 2D structure**

*TM* modes (or vice-versa).

**3.1 The general case; 3D structure**

also *Z<sup>L</sup>* as

tan *β*� = (*Z*�

If a 2D case simplifies further to become a 1D case, *Z<sup>σ</sup>* = *Z<sup>χ</sup>* = *Z* say, and the tensor for all rotations has the form [0, *Z*; −*Z*, 0]. The length *C* of the radial arm vanishes, and the diagram reduces to a pair of points on the horizontal axes, as shown in Fig. 2c.

#### **4. Invariants of rotation of the measuring axes**

It can be seen from Fig. 2a that the MT observations can be expressed in terms of seven invariants of rotation. From different sets which are possible, and evident from formal analysis (Szarka & Menvielle, 1997; Weaver et al., 2000), this chapter adopts the invariants shown in Fig. 3. Thus the eight values of the MT tensor, all of which generally change upon axes rotation, become seven invariants plus one angle. That angle, *θ*� in Fig. 2a, is the angle which defines the direction of the measuring axes (for example, with respect to north).

However a related angle for the seventh invariant, derived by Weaver et al. (2000), has the great utility that it appears again in the Mohr diagram for the phase tensor, to be introduced in Section 8 below. There, as angle *γ*, it has a simple significance concerning geologic strike.

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 87

Weaver et al. (2000) take sines of the angles to give measures in the range 0 to 1, and this technique is adopted by Marti et al. (2005). However, if ambiguity is to be avoided, such a procedure restricts the angles to being not greater than 90◦. Because examples occur for which the values of *μp*,*<sup>q</sup>* are greater than 90◦, it may be preferable to quote the invariants as angles, allowing them a 360◦ range (perhaps most usefully expressed in the range ±180◦).

**5. Some basic techniques of linear algebra applied to the analysis of the MT tensor** To gain familiarity with the MT tensor, it is instructive to explore some common steps taken

+

The second and antisymmetric part, **Za**, is immediately recognised as being of the ideal 1D

In much MT interpretation, it has been common practice to obtain "1D estimates" from 2D or 3D data by taking the average of observed *Zxy* and (−)*Zyx* as (*Zxy* − *Zyx*)/2. Equation 18 demonstrates the approximations made in such a procedure. For in taking such an average, it can be seen that the information in the first matrix term **Z<sup>s</sup>** is ignored, perhaps without justification. Diagrammatically, the procedure is equivalent, in Fig. 4a, to representing the circle on the left-hand side by the sum of the two circles on the right-hand side. The first circle on the right-hand side is then ignored, leaving the second circle which, reduced to its central

The exercise in Section 5.1 may be regarded as a separation into symmetric and 1D components, and suggests a similar separation of the observed tensor into two parts of which the second, **Z2D**, is chosen to be of ideal 2D form. The first part, **Zs2**, is found to again be

**Z** = **Z<sup>s</sup>** + **Z<sup>a</sup>** (17)

**Z** = **Zs2** + **Z2D** (19)

 <sup>0</sup> (*Zxy* <sup>−</sup> *Zyx*)/2 (*Zyx* <sup>−</sup> *Zxy*)/2 0

(18)

in matrix analysis. The steps described form part of the history of MT.

The MT tensor **Z**, if split into symmetric **Z<sup>s</sup>** and antisymmetric **Z<sup>a</sup>** parts

 *Zxx* (*Zxy* + *Zyx*)/2 (*Zxy* <sup>+</sup> *Zyx*)/2 *Zyy*

**5.1 Separation into symmetric and antisymmetric components**

**4.5 Angles or their sines?**

may be written

*Zxx Zxy Zyx Zyy*

point, is a 1D case.

symmetric.

=

MT form described in Section 3.3.

Thus the MT tensor is expressed

**5.2 Separation into symmetric and 2D components**

Fig. 3. Seven invariants of rotation for a general 3D MT tensor.

#### **4.1 Two invariants summarising the 1D character of the tensor**

The two invariants *Z<sup>L</sup> <sup>p</sup>*,*<sup>q</sup>* summarising the 1D character are straightforward and were termed "central impedances" by Lilley (1993). Note that the second could be expressed normalised by the first, and so as a ratio (and then as an angle, taking the arctangent of that ratio).

#### **4.2 Two invariants summarising the 2D character of the tensor**

Two invariants *λp*,*<sup>q</sup>* measure the 2D character, and are also straightforward. They are naturally angles, and were termed anisotropy angles by Lilley (1993). They may be expressed

$$\lambda\_{p,q} = \arcsin(\mathbb{C}\_{p,q}/Z\_{p,q}^L) \tag{16}$$

where *λp*,*<sup>q</sup>* is in the range 0 to 90◦. (Note this definition fails if *Cp*,*<sup>q</sup>* > *Z<sup>L</sup> <sup>p</sup>*,*<sup>q</sup>* and the relevant circle encloses the origin.)

#### **4.3 Two (of three) invariants summarising the 3D character of the tensor**

Two angles *μp*,*<sup>q</sup>* characterising the 3D nature of the impedance tensor are also straightforward, and are shown in Fig. 3. It may be effective to express them as their mean and difference values, because certain mechanisms for causing 3D effects, especially static distortion, give the same *μ* contribution to both the in-phase and quadrature parts of a tensor (Lilley, 1993). In such cases, static distortion of a regional 2D structure is then measured by (*μ<sup>q</sup>* + *μp*)/2, and the difference (*μ<sup>q</sup>* − *μp*) would be zero. Thus (*μ<sup>q</sup>* − *μp*), when non-zero, may be a measure of any 3D effects present, beyond static distortion.

#### **4.4 The third 3D invariant**

A third 3D invariant (*δβ* = *β<sup>q</sup>* − *βp*) can be seen from Fig. 3 to be the angle by which the two radial arms of the in-phase and quadrature circles are not parallel. It is significant as it alone, of the invariants chosen, links the in-phase and quadrature parts of an observed tensor. Another possibility for this seventh invariant is (*δβ* − *μ<sup>q</sup>* + *μp*), where the difference (*μ<sup>q</sup>* − *μp*) is removed first from *δβ*, and the departure of the radial arms from being parallel is then judged afresh.

However a related angle for the seventh invariant, derived by Weaver et al. (2000), has the great utility that it appears again in the Mohr diagram for the phase tensor, to be introduced in Section 8 below. There, as angle *γ*, it has a simple significance concerning geologic strike.

## **4.5 Angles or their sines?**

6 Will-be-set-by-IN-TECH

Z'xxq

In-phase Quadrature

L

*<sup>p</sup>*,*<sup>q</sup>* summarising the 1D character are straightforward and were termed

δβ

*<sup>p</sup>*,*q*) (16)

*<sup>p</sup>*,*<sup>q</sup>* and the relevant

Z'xyq

Z'xyp

"central impedances" by Lilley (1993). Note that the second could be expressed normalised by the first, and so as a ratio (and then as an angle, taking the arctangent of that ratio).

Two invariants *λp*,*<sup>q</sup>* measure the 2D character, and are also straightforward. They are naturally

*λp*,*<sup>q</sup>* = arcsin(*Cp*,*q*/*Z<sup>L</sup>*

Two angles *μp*,*<sup>q</sup>* characterising the 3D nature of the impedance tensor are also straightforward, and are shown in Fig. 3. It may be effective to express them as their mean and difference values, because certain mechanisms for causing 3D effects, especially static distortion, give the same *μ* contribution to both the in-phase and quadrature parts of a tensor (Lilley, 1993). In such cases, static distortion of a regional 2D structure is then measured by (*μ<sup>q</sup>* + *μp*)/2, and the difference (*μ<sup>q</sup>* − *μp*) would be zero. Thus (*μ<sup>q</sup>* − *μp*), when non-zero, may be a measure of

A third 3D invariant (*δβ* = *β<sup>q</sup>* − *βp*) can be seen from Fig. 3 to be the angle by which the two radial arms of the in-phase and quadrature circles are not parallel. It is significant as it alone, of the invariants chosen, links the in-phase and quadrature parts of an observed tensor. Another possibility for this seventh invariant is (*δβ* − *μ<sup>q</sup>* + *μp*), where the difference (*μ<sup>q</sup>* − *μp*) is removed first from *δβ*, and the departure of the radial arms from being parallel is then judged

angles, and were termed anisotropy angles by Lilley (1993). They may be expressed

where *λp*,*<sup>q</sup>* is in the range 0 to 90◦. (Note this definition fails if *Cp*,*<sup>q</sup>* > *Z<sup>L</sup>*

**4.3 Two (of three) invariants summarising the 3D character of the tensor**

μ<sup>p</sup> μ<sup>q</sup>

<sup>λ</sup>p λ<sup>q</sup>

0 0

marks the observed points

Fig. 3. Seven invariants of rotation for a general 3D MT tensor.

**4.1 Two invariants summarising the 1D character of the tensor**

**4.2 Two invariants summarising the 2D character of the tensor**

<sup>L</sup> Zq

Z'xxp

The two invariants *Z<sup>L</sup>*

circle encloses the origin.)

**4.4 The third 3D invariant**

afresh.

any 3D effects present, beyond static distortion.

Zp

Weaver et al. (2000) take sines of the angles to give measures in the range 0 to 1, and this technique is adopted by Marti et al. (2005). However, if ambiguity is to be avoided, such a procedure restricts the angles to being not greater than 90◦. Because examples occur for which the values of *μp*,*<sup>q</sup>* are greater than 90◦, it may be preferable to quote the invariants as angles, allowing them a 360◦ range (perhaps most usefully expressed in the range ±180◦).

## **5. Some basic techniques of linear algebra applied to the analysis of the MT tensor**

To gain familiarity with the MT tensor, it is instructive to explore some common steps taken in matrix analysis. The steps described form part of the history of MT.

#### **5.1 Separation into symmetric and antisymmetric components**

The MT tensor **Z**, if split into symmetric **Z<sup>s</sup>** and antisymmetric **Z<sup>a</sup>** parts

$$\mathbf{Z} = \mathbf{Z}^{\mathbf{s}} + \mathbf{Z}^{\mathbf{a}} \tag{17}$$

may be written

$$
\begin{bmatrix} Z\_{\rm xx} \ Z\_{\rm xy} \\ Z\_{\rm yx} \ Z\_{\rm yy} \end{bmatrix} = \begin{bmatrix} Z\_{\rm xx} & (Z\_{\rm xy} + Z\_{\rm yx})/2 \\ (Z\_{\rm xy} + Z\_{\rm yx})/2 & Z\_{\rm yy} \end{bmatrix} + \begin{bmatrix} 0 & (Z\_{\rm xy} - Z\_{\rm yx})/2 \\ (Z\_{\rm yx} - Z\_{\rm xy})/2 & 0 \end{bmatrix} \tag{18}
$$

The second and antisymmetric part, **Za**, is immediately recognised as being of the ideal 1D MT form described in Section 3.3.

In much MT interpretation, it has been common practice to obtain "1D estimates" from 2D or 3D data by taking the average of observed *Zxy* and (−)*Zyx* as (*Zxy* − *Zyx*)/2. Equation 18 demonstrates the approximations made in such a procedure. For in taking such an average, it can be seen that the information in the first matrix term **Z<sup>s</sup>** is ignored, perhaps without justification. Diagrammatically, the procedure is equivalent, in Fig. 4a, to representing the circle on the left-hand side by the sum of the two circles on the right-hand side. The first circle on the right-hand side is then ignored, leaving the second circle which, reduced to its central point, is a 1D case.

#### **5.2 Separation into symmetric and 2D components**

The exercise in Section 5.1 may be regarded as a separation into symmetric and 1D components, and suggests a similar separation of the observed tensor into two parts of which the second, **Z2D**, is chosen to be of ideal 2D form. The first part, **Zs2**, is found to again be symmetric.

Thus the MT tensor is expressed

$$\mathbf{Z} = \mathbf{Z}^{\mathbf{s}2} + \mathbf{Z}^{\mathbf{2D}} \tag{19}$$

**5.3.1 Eigenvalues are conjugate pairs**

**5.3.2 Eigenvalues are real and equal**

**5.3.3 Eigenvalues are real and different**

The product of the two eigenvalues is given by

**5.4 Eigenvalues on Mohr diagrams**

and is positive if det**Z** is positive, and negative if det**Z** is negative.

*xy* = 0 and *Z*�

The third case occurs when

respectively of **Z**.

of the axes.

them do not exist.

both vertical axes, *Z*�

Again, the product of the two roots, *ζ*1*ζ*2, will be positive.

The second case occurs when

In fact

(*Zxx* <sup>+</sup> *Zyy*)<sup>2</sup> <sup>+</sup> <sup>4</sup>(*ZxyZyx* <sup>−</sup> *ZxxZyy*) <sup>&</sup>lt; <sup>0</sup> (23)

(*Zxx* <sup>+</sup> *Zyy*)<sup>2</sup> <sup>+</sup> <sup>4</sup>(*ZxyZyx* <sup>−</sup> *ZxxZyy*) = <sup>0</sup> (25)

(*Zxx* <sup>+</sup> *Zyy*)<sup>2</sup> <sup>+</sup> <sup>4</sup>(*ZxyZyx* <sup>−</sup> *ZxxZyy*) <sup>&</sup>gt; <sup>0</sup> (27)

*ζ*1*ζ*<sup>2</sup> = det **Z** (28)

*ζ*<sup>1</sup> = *ζ*<sup>2</sup> = (*Zxx* + *Zyy*)/2 (26)

(*Zxx* <sup>−</sup> *Zyy*)<sup>2</sup> <sup>+</sup> <sup>4</sup>*ZxyZyx* <sup>&</sup>lt; <sup>0</sup> (24)

and the two roots of the conjugate equation form a conjugate pair. The product of the two roots, *ζ*1*ζ*2, will always be positive. An equivalent way of expressing Inequality 23 is as

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 89

The two roots of the conjugate equation are now both real (positive or negative), and equal.

The two roots of the conjugate equation are now both real, different, and positive or negative depending on the signs of (*Zxx* + *Zyy*) and (*ZxxZyy* − *ZxyZyx*), the trace and determinant

The three eigenvalue cases discussed in the preceding three subsections are clearly defined when the data are plotted on Mohr diagrams. Eigenvalues which are real are shown graphically. The directions of their eigenvectors may be read from the diagrams remembering, in Fig. 2a, the 2*θ*� anticlockwise rotation of the radial arm for, in Fig. 1, the *θ*� clockwise rotation

Thus for the case of Section 5.3.1, circles which (as in Fig. 5a) do not touch the vertical axes obey Inequality 23. Their eigenvalues are complex conjugate pairs, and real eigenvectors for

Secondly, for the case discussed in Section 5.3.2, Fig. 5b shows a circle which is just touching

*yx* = 0. The eigenvalues may be read off the *Z*�

*xx* axis, as

The first case occurs when

Fig. 4. Separation ofa2x2 matrix, a: into symmetric and 1D parts as in Equation 18; b: into symmetric and 2D parts as in Equation 20. The diagrams are drawn for **Z** = [1.75, 3.30; −0.70, 0.25]. Thus case a is: **Z** = [1.75, 1.3; 1.3, 0.25]+[0, 2.0; −2.0, 0], and case b is: **Z** = [1.0, 0; 0, 1.0]+[0.75, 3.3; −0.7, −0.75].

and by partition as

$$
\begin{bmatrix} Z\_{\text{xx}} \ Z\_{\text{xy}} \\ Z\_{\text{yx}} \ Z\_{\text{yy}} \end{bmatrix} = \begin{bmatrix} (\mathbf{Z}\_{\text{xx}} + \mathbf{Z}\_{\text{yy}})/2 & \mathbf{0} \\ \mathbf{0} & (\mathbf{Z}\_{\text{xx}} + \mathbf{Z}\_{\text{yy}})/2 \end{bmatrix} + \begin{bmatrix} (\mathbf{Z}\_{\text{xx}} - \mathbf{Z}\_{\text{yy}})/2 & \mathbf{Z}\_{\text{xy}} \\ \mathbf{Z}\_{\text{yx}} & (-\mathbf{Z}\_{\text{xx}} + \mathbf{Z}\_{\text{yy}})/2 \end{bmatrix} \tag{20}
$$

where the second part is now of ideal 2D form.

Equation 20 is expressed in Mohr diagrams in Fig. 4b, where it can be seen that the first part is a point on the vertical axis, and the second part, **Z2D**, taken by itself plots as an ideal 2D circle with centre on the horizontal *Z*� *xy* axis. The intersections of the circle with the axis give the values of the *TE* and *TM* impedances for this (now artificially) ideal tensor.

Thus the common practice with MT data, when seeking *TE* and *TM* values, of finding the maximum and minimum values of *Z*� *xy* as the observing axes are rotated, can be seen to be tantamount to an assumption that the **Zs2** part of the matrix be ignored, so that the circle for the **Z2D** part indeed plots in ideal 2D form.

#### **5.3 Eigenvalue analysis**

Eigenvalues *ζ*<sup>1</sup> and *ζ*<sup>2</sup> of a matrix **Z** are found by solving the characteristic equation

$$\left(\zeta^2 - (\mathbf{Z}\_{xx} + \mathbf{Z}\_{yy})\zeta + \mathbf{Z}\_{xx}\mathbf{Z}\_{yy} - \mathbf{Z}\_{xy}\mathbf{Z}\_{yx} = \mathbf{0}\tag{21}$$

to obtain

$$\mathcal{L}\_1 \mathcal{L}\_2 = (\mathcal{Z}\_{\text{xx}} + \mathcal{Z}\_{yy})/2 \pm \left[ (\mathcal{Z}\_{\text{xx}} + \mathcal{Z}\_{yy})^2 + 4(\mathcal{Z}\_{\text{xy}}\mathcal{Z}\_{yx} - \mathcal{Z}\_{\text{xx}}\mathcal{Z}\_{yy}) \right]^\frac{1}{2}/2 \tag{22}$$

Three cases are possible and of interest. Each case will be discussed separately.

#### **5.3.1 Eigenvalues are conjugate pairs**

The first case occurs when

8 Will-be-set-by-IN-TECH

Fig. 4. Separation ofa2x2 matrix, a: into symmetric and 1D parts as in Equation 18; b: into

**Z** = [1.75, 3.30; −0.70, 0.25]. Thus case a is: **Z** = [1.75, 1.3; 1.3, 0.25]+[0, 2.0; −2.0, 0], and case

 + 

Equation 20 is expressed in Mohr diagrams in Fig. 4b, where it can be seen that the first part is a point on the vertical axis, and the second part, **Z2D**, taken by itself plots as an ideal 2D

Thus the common practice with MT data, when seeking *TE* and *TM* values, of finding the

tantamount to an assumption that the **Zs2** part of the matrix be ignored, so that the circle for

+

(*Zxx* − *Zyy*)/2 *Zxy*

*xy* axis. The intersections of the circle with the axis give

*xy* as the observing axes are rotated, can be seen to be

*<sup>ζ</sup>*<sup>2</sup> <sup>−</sup> (*Zxx* <sup>+</sup> *Zyy*)*<sup>ζ</sup>* <sup>+</sup> *ZxxZyy* <sup>−</sup> *ZxyZyx* <sup>=</sup> <sup>0</sup> (21)

*Zyx* (−*Zxx* + *Zyy*)/2

<sup>2</sup> /2 (22)

(20)

+

=

symmetric and 2D parts as in Equation 20. The diagrams are drawn for

0 (*Zxx* + *Zyy*)/2

the values of the *TE* and *TM* impedances for this (now artificially) ideal tensor.

Eigenvalues *ζ*<sup>1</sup> and *ζ*<sup>2</sup> of a matrix **Z** are found by solving the characteristic equation

Three cases are possible and of interest. Each case will be discussed separately.

*<sup>ζ</sup>*1, *<sup>ζ</sup>*<sup>2</sup> = (*Zxx* <sup>+</sup> *Zyy*)/2 <sup>±</sup> [(*Zxx* <sup>+</sup> *Zyy*)<sup>2</sup> <sup>+</sup> <sup>4</sup>(*ZxyZyx* <sup>−</sup> *ZxxZyy*)] <sup>1</sup>

b is: **Z** = [1.0, 0; 0, 1.0]+[0.75, 3.3; −0.7, −0.75].

where the second part is now of ideal 2D form.

circle with centre on the horizontal *Z*�

maximum and minimum values of *Z*�

the **Z2D** part indeed plots in ideal 2D form.

(*Zxx* + *Zyy*)/2 0

=

a

b

and by partition as

 = 

**5.3 Eigenvalue analysis**

to obtain

*Zxx Zxy Zyx Zyy*

$$((Z\_{xx} + Z\_{yy})^2 + 4(Z\_{xy}Z\_{yx} - Z\_{xx}Z\_{yy}) < 0\tag{23}$$

and the two roots of the conjugate equation form a conjugate pair. The product of the two roots, *ζ*1*ζ*2, will always be positive. An equivalent way of expressing Inequality 23 is as

$$\left(\left(Z\_{xx} - Z\_{yy}\right)^2 + 4Z\_{xy}Z\_{yx} < 0\right) \tag{24}$$

#### **5.3.2 Eigenvalues are real and equal**

The second case occurs when

$$(Z\_{xx} + Z\_{yy})^2 + 4(Z\_{xy}Z\_{yx} - Z\_{xx}Z\_{yy}) = 0\tag{25}$$

The two roots of the conjugate equation are now both real (positive or negative), and equal. In fact

$$\mathbb{Z}\_1 = \mathbb{Z}\_2 = (\mathbb{Z}\_{xx} + \mathbb{Z}\_{yy})/2 \tag{26}$$

Again, the product of the two roots, *ζ*1*ζ*2, will be positive.

#### **5.3.3 Eigenvalues are real and different**

The third case occurs when

$$((Z\_{xx} + Z\_{yy})^2 + 4(Z\_{xy}Z\_{yx} - Z\_{xx}Z\_{yy}) > 0\tag{27}$$

The two roots of the conjugate equation are now both real, different, and positive or negative depending on the signs of (*Zxx* + *Zyy*) and (*ZxxZyy* − *ZxyZyx*), the trace and determinant respectively of **Z**.

The product of the two eigenvalues is given by

$$
\mathbb{Z}\_1 \mathbb{Z}\_2 = \det \mathbf{Z} \tag{28}
$$

and is positive if det**Z** is positive, and negative if det**Z** is negative.

#### **5.4 Eigenvalues on Mohr diagrams**

The three eigenvalue cases discussed in the preceding three subsections are clearly defined when the data are plotted on Mohr diagrams. Eigenvalues which are real are shown graphically. The directions of their eigenvectors may be read from the diagrams remembering, in Fig. 2a, the 2*θ*� anticlockwise rotation of the radial arm for, in Fig. 1, the *θ*� clockwise rotation of the axes.

Thus for the case of Section 5.3.1, circles which (as in Fig. 5a) do not touch the vertical axes obey Inequality 23. Their eigenvalues are complex conjugate pairs, and real eigenvectors for them do not exist.

Secondly, for the case discussed in Section 5.3.2, Fig. 5b shows a circle which is just touching both vertical axes, *Z*� *xy* = 0 and *Z*� *yx* = 0. The eigenvalues may be read off the *Z*� *xx* axis, as

**6. Singular value decomposition**

in-phase and quadrature parts as

2D tensor, for example by expressing Equation 29 as

*<sup>θ</sup>ep* <sup>=</sup> <sup>1</sup> 2 arctan

*<sup>θ</sup>hp* <sup>=</sup> <sup>1</sup> 2 arctan

**E***p*,*<sup>q</sup>* = **R**(−*θep*,*q*)

*Zxx <sup>p</sup> Zxy <sup>p</sup> Zyx <sup>p</sup>*

*Zyy <sup>p</sup>*

of the tensor separately.

**W**.

is factored into

where

and

*q* replacing *p*.

It is instructive also to apply SVD to an MT tensor, taking the in-phase and quadrature parts

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 91

where the columns of **V** are eigenvectors of **ATA**, and **A<sup>T</sup>** denotes the transpose of **A**. The columns of **U** (which are eigenvectors of **AAT**) may be found by multiplying **A** by the columns of **V**. The singular values on the diagonal of **Σ** are the square roots (taken positive by convention, a most important point in the present context) of the non-zero eigenvalues of **AAT**. As **Σ** is diagonal, it is straightforward to convert it to the antidiagonal form of an ideal

where **<sup>W</sup>** = [0, 1; <sup>−</sup>1, 0] and **<sup>V</sup><sup>T</sup>** is pre-multiplied by **<sup>W</sup>**−**1**, with **<sup>W</sup>**−**<sup>1</sup>** denoting the inverse of

Equations for this form of the SVD of an MT tensor were derived in an earlier paper (Lilley, 1998). Taking phase relative to the **H** signal, so that **H***<sup>q</sup>* is zero, Equation 1 is written in its

where the electric and magnetic observation axes are rotated clockwise independently, the electric axes by *θep*,*<sup>q</sup>* and the magnetic axes by *θhp*,*q*. Thus the in-phase part of the tensor, **Z***p*,

= **R**(−*θep*)

*Zyy <sup>p</sup>* − *Zxx <sup>p</sup> Zxy <sup>p</sup>* + *Zyx <sup>p</sup>*

*Zyy <sup>p</sup>* − *Zxx <sup>p</sup> Zxy <sup>p</sup>* + *Zyx <sup>p</sup>*

 0 Υ*p*,*<sup>q</sup>* −Ψ*p*,*<sup>q</sup>* 0

 0 Υ*<sup>p</sup>* −Ψ*<sup>p</sup>* 0

+ arctan

− arctan

Υ*<sup>p</sup>* − Ψ*<sup>p</sup>* = cos(*θep* + *θhp*)[(*Zxy <sup>p</sup>* + *Zyx <sup>p</sup>*) − tan(*θep* + *θhp*)(*Zxx <sup>p</sup>* − *Zyy <sup>p</sup>*)] (35)

Υ*<sup>p</sup>* + Ψ*<sup>p</sup>* = cos(*θep* − *θhp*)[(*Zxy <sup>p</sup>* − *Zyx <sup>p</sup>*) + tan(*θep* − *θhp*)(*Zxx <sup>p</sup>* + *Zyy <sup>p</sup>*)] (36)

(Equation 36 corrects equation 22 of Lilley (1998), where a negative sign is missing.) A similar set to Equations 32, 33, 34, 35 and 36 applies for the quadrature part of a tensor, with subscript

The quantities Υ*p*,*<sup>q</sup>* and Ψ*p*,*<sup>q</sup>* are termed principal values. An examination of the possibility

that one or both of Ψ*p*,*<sup>q</sup>* may be zero or negative is addressed in Section 13.3.

*Zyy <sup>p</sup>* + *Zxx <sup>p</sup> Zxy <sup>p</sup>* − *Zyx <sup>p</sup>*

*Zyy <sup>p</sup>* + *Zxx <sup>p</sup> Zxy <sup>p</sup>* − *Zyx <sup>p</sup>*

**A** = **UΣV<sup>T</sup>** (29)

**A** = **UΣWW**<sup>−</sup>**1VT** (30)

**R**(*θhp*,*q*)**H** (31)

**R**(*θhp*) (32)

(33)

(34)

As described for example by Strang (2005), decomposition by SVD factors a matrix **A** into

Fig. 5. Eigenvector directions (according to the radial arrows) and eigenvalues (where the arrowheads touch an axis) displayed on Mohr diagrams fora2x2 matrix. a: Eigenvalues are complex conjugates, and are not evident on the diagram. b: Eigenvalues are real and equal; both eigenvectors are the same. c: Eigenvalues are real and different; eigenvectors are not orthogonal. d: Eigenvalues are real and different; eigenvectors are orthogonal. e: Associated MT eigenvalues are real and different; eigenvectors are not orthogonal. f: Associated MT eigenvalues are real and different; eigenvectors are orthogonal (the 2D case). In Fig. 5a axes are also drawn for *Z yx* and *Z yy*. Between the *Z xx* and *Z yy* axes, the product *Z xy*.*Z yx* < 0. On the *Z xx* and *Z yy* axes, *Z xy*.*Z yx* = 0. Outside the *Z xx* and *Z yy* axes, *Z xy*.*Z yx* > 0.

(*Zxx* + *Zyy*)/2. The direction of the repeated eigenvector corresponds to the direction of a radial arm which is horizontal in the diagram, as shown.

Thirdly, for the case discussed in Section 5.3.3, in Fig. 5c a circle is shown which, like examples to be discussed below, now crosses the vertical axes. The two eigenvalues may again be read off the *Z xx* axis where the circle cuts this axis, and the two eigenvector directions are given by the *θ* values for the radial arms to these points.

The example in Fig. 5c demonstrates that real eigenvalues correspond to the MT case of "phase going out of quadrant". Also it can be seen, by an extension of the discussion in Section 5.3.3 above, that for a Mohr circle to not capture the origin the product of the two eigenvalues must be positive; i.e. det**Z** must be positive.

Of equal interest in MT are the eigenvectors of the associated problem, which correspond to intersections of a circle with the horizontal *Z xy* axis. As shown in Fig. 5e, these might be regarded as an approximation to 2D *TE* and *TM* values, and for the true 2D case (Fig. 5f) they will indeed be so. These eigenvalues may be found by similar formalism, or by elementary trigonometry based on Fig. 2a. They may be evident from inspection, as in Lilley (1993).

#### **6. Singular value decomposition**

10 Will-be-set-by-IN-TECH

a b c

d e f

*yy*. Between the *Z*

*yx* = 0. Outside the *Z*

are also drawn for *Z*

*xx* and *Z*

the *Z*

off the *Z*

*yx* and *Z*

the *θ* values for the radial arms to these points.

intersections of a circle with the horizontal *Z*

*xy*.*Z*

radial arm which is horizontal in the diagram, as shown.

*yy* axes, *Z*

be positive; i.e. det**Z** must be positive.

Fig. 5. Eigenvector directions (according to the radial arrows) and eigenvalues (where the arrowheads touch an axis) displayed on Mohr diagrams fora2x2 matrix. a: Eigenvalues are complex conjugates, and are not evident on the diagram. b: Eigenvalues are real and equal; both eigenvectors are the same. c: Eigenvalues are real and different; eigenvectors are not orthogonal. d: Eigenvalues are real and different; eigenvectors are orthogonal. e: Associated MT eigenvalues are real and different; eigenvectors are not orthogonal. f: Associated MT eigenvalues are real and different; eigenvectors are orthogonal (the 2D case). In Fig. 5a axes

*xx* and *Z*

*xx* axis where the circle cuts this axis, and the two eigenvector directions are given by

(*Zxx* + *Zyy*)/2. The direction of the repeated eigenvector corresponds to the direction of a

Thirdly, for the case discussed in Section 5.3.3, in Fig. 5c a circle is shown which, like examples to be discussed below, now crosses the vertical axes. The two eigenvalues may again be read

The example in Fig. 5c demonstrates that real eigenvalues correspond to the MT case of "phase going out of quadrant". Also it can be seen, by an extension of the discussion in Section 5.3.3 above, that for a Mohr circle to not capture the origin the product of the two eigenvalues must

Of equal interest in MT are the eigenvectors of the associated problem, which correspond to

regarded as an approximation to 2D *TE* and *TM* values, and for the true 2D case (Fig. 5f) they will indeed be so. These eigenvalues may be found by similar formalism, or by elementary trigonometry based on Fig. 2a. They may be evident from inspection, as in Lilley (1993).

*xx* and *Z*

*yy* axes, the product *Z*

*xy*.*Z*

*xy* axis. As shown in Fig. 5e, these might be

*yx* > 0.

*yy* axes, *Z*

*xy*.*Z*

*yx* < 0. On

It is instructive also to apply SVD to an MT tensor, taking the in-phase and quadrature parts of the tensor separately.

As described for example by Strang (2005), decomposition by SVD factors a matrix **A** into

$$\mathbf{A} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{\mathrm{T}} \tag{29}$$

where the columns of **V** are eigenvectors of **ATA**, and **A<sup>T</sup>** denotes the transpose of **A**. The columns of **U** (which are eigenvectors of **AAT**) may be found by multiplying **A** by the columns of **V**. The singular values on the diagonal of **Σ** are the square roots (taken positive by convention, a most important point in the present context) of the non-zero eigenvalues of **AAT**. As **Σ** is diagonal, it is straightforward to convert it to the antidiagonal form of an ideal 2D tensor, for example by expressing Equation 29 as

$$\mathbf{A} = \mathbf{U}\Sigma\mathbf{W}\mathbf{W}^{-1}\mathbf{V}^{T} \tag{30}$$

where **<sup>W</sup>** = [0, 1; <sup>−</sup>1, 0] and **<sup>V</sup><sup>T</sup>** is pre-multiplied by **<sup>W</sup>**−**1**, with **<sup>W</sup>**−**<sup>1</sup>** denoting the inverse of **W**.

Equations for this form of the SVD of an MT tensor were derived in an earlier paper (Lilley, 1998). Taking phase relative to the **H** signal, so that **H***<sup>q</sup>* is zero, Equation 1 is written in its in-phase and quadrature parts as

$$\mathbf{E}\_{p,q} = \mathbf{R}(-\theta e\_{p,q}) \begin{bmatrix} 0 & \mathbf{Y}\_{p,q} \\ -\mathbf{Y}\_{p,q} & 0 \end{bmatrix} \mathbf{R}(\theta h\_{p,q}) \mathbf{H} \tag{31}$$

where the electric and magnetic observation axes are rotated clockwise independently, the electric axes by *θep*,*<sup>q</sup>* and the magnetic axes by *θhp*,*q*. Thus the in-phase part of the tensor, **Z***p*, is factored into

$$
\begin{bmatrix} Z\_{\mathbf{x}\mathbf{x}p} Z\_{\mathbf{x}y\_p} \\ Z\_{\mathbf{y}\mathbf{x}\_p} Z\_{\mathbf{y}y\_p} \end{bmatrix} = \mathbf{R}(-\theta e\_p) \begin{bmatrix} 0 & \mathbf{Y}\_p \\ -\mathbf{Y}\_p & 0 \end{bmatrix} \mathbf{R}(\theta \mathbf{h}\_p) \tag{32}
$$

where

$$\theta e\_p = \frac{1}{2} \left[ \arctan \frac{Z\_{yy\_p} - Z\_{xx\_p}}{Z\_{xy\_p} + Z\_{yx\_p}} + \arctan \frac{Z\_{yy\_p} + Z\_{xx\_p}}{Z\_{xy\_p} - Z\_{yx\_p}} \right] \tag{33}$$

$$\theta h\_p = \frac{1}{2} \left[ \arctan \frac{Z\_{yy\_p} - Z\_{xx\_p}}{Z\_{xy\_p} + Z\_{yx\_p}} - \arctan \frac{Z\_{yy\_p} + Z\_{xx\_p}}{Z\_{xy\_p} - Z\_{yx\_p}} \right] \tag{34}$$

$$(\Psi\_p - \Psi\_p = \cos(\theta e\_p + \theta h\_p)[(Z\_{xy\_p} + Z\_{yx\_p}) - \tan(\theta e\_p + \theta h\_p)(Z\_{xx\_p} - Z\_{yy\_p})] \tag{35}$$

and

$$\Psi\_p + \Psi\_p = \cos(\theta e\_p - \theta h\_p)[(Z\_{xy} - Z\_{yx\_p}) + \tan(\theta e\_p - \theta h\_p)(Z\_{xxp} + Z\_{yy\_p})] \tag{36}$$

(Equation 36 corrects equation 22 of Lilley (1998), where a negative sign is missing.) A similar set to Equations 32, 33, 34, 35 and 36 applies for the quadrature part of a tensor, with subscript *q* replacing *p*.

The quantities Υ*p*,*<sup>q</sup>* and Ψ*p*,*<sup>q</sup>* are termed principal values. An examination of the possibility that one or both of Ψ*p*,*<sup>q</sup>* may be zero or negative is addressed in Section 13.3.

and then as

**Z***<sup>p</sup>* = 

the values of Υ*<sup>p</sup>* and Ψ*p*.

field changes.

cos 31.71◦ − sin 31.71◦ sin 31.71◦ cos 31.71◦

and *u*<sup>2</sup> which are as given for *E*�

possible if regional anisotropy is high.

there is some rotation of axes for which both *Z*�

0 8.0902

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 93

The columns of the first matrix [cos 31.71◦, − sin 31.71◦; sin 31.71◦, cos 31.71◦], which have been derived from **U** in Equation 29, can be seen to define unit vectors *u*<sup>1</sup>

matrix [cos 21.41◦, sin 21.41◦; <sup>−</sup> sin 21.41◦, cos 21.41◦], which have been derived from **<sup>V</sup><sup>T</sup>** in

Also, as is evident, the singular values in the diagonal of matrix **Σ** (the second matrix) give

Note that by rotation of the electric and magnetic axes separately, the MT tensor has been reduced to an ideal 2D form. In a search for the nearest 2D model in a 3D situation the results of SVD may be valuable to bear in mind; the magnetic axes may be an indication of regional strike, with the electric axes showing the brunt of the distortion. A disadvantage of such an analysis however is that a rather artificial conductivity structure is required to simply twist the electric field at an observing site to explain the different rotations required of the *E* and *H* axes. A more specific (if simple) model for local distortion is widely accepted, and will be discussed in Section 11. With this model, in the case of near-singular MT data, it will be shown that a surficial strike direction is determined. However regional strike determination remains

It is common experience to find very strong anisotropy in an observed tensor, both for distorted 2D cases, and indeed generally. As a consequence, the tensor approaches a condition of singularity, in both its in-phase and quadrature parts. In a Mohr diagram the condition of singularity is shown by a circle touching the origin. If for example **Z***<sup>p</sup>* is a singular tensor, then

zero, when error is taken into account). For the student of linear algebra, an example of a null space occurs: it is the line of the direction of nil electric field change, holding for all magnetic

A condition number may be used to warn that singularity is being approached (Strang, 2005). When the condition number becomes high in some sense, the matrix is said to be ill-conditioned (Press et al., 1989). The condition number suggested by Strang (2005) is the norm of the matrix (sometimes called the spectral norm) multiplied by the norm of the inverse of the matrix; or equivalently, the greater principal value of the matrix divided by the lesser

The greater and lesser principal values Υ*p*,*<sup>q</sup>* and Ψ*p*,*<sup>q</sup>* are given by Equations 35 and 36 and, as discussed above, are the singular values of the matrix. Following the convention that singular values are never negative, a modulus sign is put into the denominator of Equation 41 to cover cases where computation of the lesser principal value produces a negative number. In terms

principal value. For the2x2 matrix **Z** in Equation 1 the condition numbers *κp*,*<sup>q</sup>* are

*xx <sup>p</sup>* and *Z*�

*<sup>x</sup>* and *E*�

Equation 29, can be seen to define unit vectors *v*<sup>1</sup> and *v*<sup>2</sup> as given for *H*�

**7. A condition number to measure singularity in an MT tensor**

<sup>−</sup>3.0902 0 cos 21.41◦ sin 21.41◦

− sin 21.41◦ cos 21.41◦

*<sup>y</sup>* in Fig. 6. The rows of the third

*xy <sup>p</sup>* are zero (or indistinguishable from

*κp*,*<sup>q</sup>* = Υ*p*,*q*/|Ψ*p*,*q*| (41)

*<sup>x</sup>* and *H*�

(40)

*<sup>y</sup>* in Fig. 6.

Fig. 6. Diagram showing the Mohr representation and the SVD of the matrix in Equation 37. The values of *θep*, *θhp*, Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* are evident as 32◦, 21◦, 8.1 and 3.1. The axes for *E*� *<sup>x</sup>*, *E*� *<sup>y</sup>* and *H*� *<sup>x</sup>*, *H*� *<sup>y</sup>* show the rotations, from *Ex*, *Ey* and *Hx*, *Hy* by *θep* and *θhp* respectively, to give an ideal 2D antidiagonal response.

The quantities arising from the SVD may be displayed on a Mohr diagram as in Fig. 6, which has been drawn for a tensor

$$\mathbf{Z}\_p = [-1, 7; -4, 3] \tag{37}$$

The values of *θep*, *θhp*, Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* may be evaluated by the equations above as 31.7◦, 21.4◦, 8.09 and 3.09 respectively. These values may also be read off the figure.

When the matrix of Equation 37 is put into a standard computing routine, the SVD returned is commonly in the form of Equation 29:

$$\mathbf{Z}\_p = \begin{bmatrix} -.8507 \ -.5257 \\ -.5257 \ .8507 \end{bmatrix} \begin{bmatrix} 8.0902 & 0 \\ 0 & 3.0902 \end{bmatrix} \begin{bmatrix} .3651 \ -.9310 \\ -.9310 \ -.3651 \end{bmatrix} \tag{38}$$

which may be given the form of Equation 32 by expressing it first as

$$\mathbf{Z}\_p = \begin{bmatrix} .8507 \ -.5257 \\ .5257 \ .8507 \end{bmatrix} \begin{bmatrix} 0 & 8.0902 \\ -3.0902 & 0 \end{bmatrix} \begin{bmatrix} .9310 & .3651 \\ -3651 & .9310 \end{bmatrix} \tag{39}$$

and then as

12 Will-be-set-by-IN-TECH

0 5.5

u1 v1

u2 v2

Fig. 6. Diagram showing the Mohr representation and the SVD of the matrix in Equation 37.

*<sup>y</sup>* show the rotations, from *Ex*, *Ey* and *Hx*, *Hy* by *θep* and *θhp* respectively, to give an

The quantities arising from the SVD may be displayed on a Mohr diagram as in Fig. 6, which

The values of *θep*, *θhp*, Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* may be evaluated by the equations above as 31.7◦, 21.4◦, 8.09

When the matrix of Equation 37 is put into a standard computing routine, the SVD returned

8.0902 0

 0 8.0902 −3.0902 0

0 3.0902

Ey

Ey'

The values of *θep*, *θhp*, Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* are evident as 32◦, 21◦, 8.1 and 3.1. The axes for *E*�

θep - θhp

Ex Ex'

θep

and 3.09 respectively. These values may also be read off the figure.

−.8507 −.5257 −.5257 .8507

.8507 −.5257 .5257 .8507

which may be given the form of Equation 32 by expressing it first as

Ψp

Υp

Z'xxp

1

*H*� *<sup>x</sup>*, *H*�

ideal 2D antidiagonal response.

is commonly in the form of Equation 29:

**Z***<sup>p</sup>* = 

> **Z***<sup>p</sup>* =

has been drawn for a tensor

observed point (Zxxp =-1, Zxyp =7)

Hy

*<sup>x</sup>*, *E*� *<sup>y</sup>* and

(38)

(39)

Hy'

**Z***<sup>p</sup>* = [−1, 7; −4, 3] (37)

 .3651 <sup>−</sup>.9310 −.9310 −.3651

 .9310 .3651 −.3651 .9310 2θep

Hx Hx'

θ hp

Z'xyp

$$\mathbf{Z}\_p = \begin{bmatrix} \cos 31.71^\circ - \sin 31.71^\circ \\ \sin 31.71^\circ \end{bmatrix} \begin{bmatrix} 0 & 8.0902 \\ -3.0902 & 0 \end{bmatrix} \begin{bmatrix} \cos 21.41^\circ & \sin 21.41^\circ \\ -\sin 21.41^\circ \cos 21.41^\circ \end{bmatrix} \tag{40}$$

The columns of the first matrix [cos 31.71◦, − sin 31.71◦; sin 31.71◦, cos 31.71◦], which have been derived from **U** in Equation 29, can be seen to define unit vectors *u*<sup>1</sup> and *u*<sup>2</sup> which are as given for *E*� *<sup>x</sup>* and *E*� *<sup>y</sup>* in Fig. 6. The rows of the third matrix [cos 21.41◦, sin 21.41◦; <sup>−</sup> sin 21.41◦, cos 21.41◦], which have been derived from **<sup>V</sup><sup>T</sup>** in Equation 29, can be seen to define unit vectors *v*<sup>1</sup> and *v*<sup>2</sup> as given for *H*� *<sup>x</sup>* and *H*� *<sup>y</sup>* in Fig. 6. Also, as is evident, the singular values in the diagonal of matrix **Σ** (the second matrix) give the values of Υ*<sup>p</sup>* and Ψ*p*.

Note that by rotation of the electric and magnetic axes separately, the MT tensor has been reduced to an ideal 2D form. In a search for the nearest 2D model in a 3D situation the results of SVD may be valuable to bear in mind; the magnetic axes may be an indication of regional strike, with the electric axes showing the brunt of the distortion. A disadvantage of such an analysis however is that a rather artificial conductivity structure is required to simply twist the electric field at an observing site to explain the different rotations required of the *E* and *H* axes. A more specific (if simple) model for local distortion is widely accepted, and will be discussed in Section 11. With this model, in the case of near-singular MT data, it will be shown that a surficial strike direction is determined. However regional strike determination remains possible if regional anisotropy is high.

#### **7. A condition number to measure singularity in an MT tensor**

It is common experience to find very strong anisotropy in an observed tensor, both for distorted 2D cases, and indeed generally. As a consequence, the tensor approaches a condition of singularity, in both its in-phase and quadrature parts. In a Mohr diagram the condition of singularity is shown by a circle touching the origin. If for example **Z***<sup>p</sup>* is a singular tensor, then there is some rotation of axes for which both *Z*� *xx <sup>p</sup>* and *Z*� *xy <sup>p</sup>* are zero (or indistinguishable from zero, when error is taken into account). For the student of linear algebra, an example of a null space occurs: it is the line of the direction of nil electric field change, holding for all magnetic

field changes.

A condition number may be used to warn that singularity is being approached (Strang, 2005). When the condition number becomes high in some sense, the matrix is said to be ill-conditioned (Press et al., 1989). The condition number suggested by Strang (2005) is the norm of the matrix (sometimes called the spectral norm) multiplied by the norm of the inverse of the matrix; or equivalently, the greater principal value of the matrix divided by the lesser principal value. For the2x2 matrix **Z** in Equation 1 the condition numbers *κp*,*<sup>q</sup>* are

$$\kappa\_{p,q} = \mathbf{Y}\_{p,q} / |\mathbf{Y}\_{p,q}| \tag{41}$$

The greater and lesser principal values Υ*p*,*<sup>q</sup>* and Ψ*p*,*<sup>q</sup>* are given by Equations 35 and 36 and, as discussed above, are the singular values of the matrix. Following the convention that singular values are never negative, a modulus sign is put into the denominator of Equation 41 to cover cases where computation of the lesser principal value produces a negative number. In terms

Fig. 7. The condition number *κp*,*<sup>q</sup>* displayed as a function of the anisotropy angle *λp*,*<sup>q</sup>* in degrees, using linear (left) and logarithmic (right) scales for *κp*,*q*.

of the Mohr representations in Figs 2, 3 and 6, from Equation 41 the condition numbers are given by

$$\kappa\_{p\mathcal{A}} = (\mathbf{Z}\_{p,q}^{\mathcal{L}} + \mathbb{C}\_{p,q}) / |(\mathbf{Z}\_{p,q}^{\mathcal{L}} - \mathbb{C}\_{p,q})| \tag{42}$$

that is

$$\kappa\_{p,q} = 2/(\csc \lambda\_{p,q} - 1) + 1 \tag{43}$$

remembering Equation 16. Fig. 7 shows *κp*,*<sup>q</sup>* as a function of *λp*,*q*. While *λp*,*<sup>q</sup>* itself is a measure of condition, Fig. 7 demonstrates that *κp*,*<sup>q</sup>* is more sensitive than *λp*,*<sup>q</sup>* as ill-condition is approached. Note that *κp*,*<sup>q</sup>* ≥ 1.

An example is given in Fig. 8 of condition numbers determined for a set of data from a recent MT site, NQ142, in north Queensland, Australia. An increase in condition number above 10 monitors the decrease of the lesser principle value Ψ*p*,*q*, which becomes negligible. The variation of angle *θep*,*<sup>q</sup>* with period (T) stabilizes as condition numbers rise, an effect discussed below in Section 11.

#### **8. The phase tensor**

Caldwell et al. (2004) introduced a "phase tensor" based on the magnetotelluric tensor, and the concept was further developed by Bibby et al. (2005). The phase tensor is a real matrix **Φ**, defined by

$$
\Phi = \mathbf{Z\_p^{-1} Z\_q} \tag{44}
$$




0

Z'xxp T1/2

50


10-4 10-3 10-2 10-1 10<sup>0</sup> 101 10<sup>2</sup> 10<sup>3</sup> 10<sup>4</sup> 10<sup>5</sup> T (s)

circumstance; however condition numbers are computed nevertheless.




0

Z'xxq T1/2

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 95

50


10-4 10-3 10-2 10-1 10<sup>0</sup> 10<sup>1</sup> 10<sup>2</sup> 10<sup>3</sup> 104 10<sup>5</sup> T (s)

κ<sup>q</sup>

Fig. 8. The NQ142 data plotted as Mohr diagrams and analysed by SVD to give: angles *θep*, *θeq*, *θhp* and *θhq*; the principal values Υ*p*, Ψ*p*, Υ*<sup>q</sup>* and Ψ*q*; and the condition numbers *κ<sup>p</sup>* and *κq*. The variation of period with colour in the circles is the same as for the other plots. The Mohr diagrams follow the form of Fig. 2a, and have their impedance values multiplied by the square root of period (T) to make the plots more compact. Where circles enclose the origin, their lesser principal value (Ψ) is given an artificially high value to flag this

Ψq ,

Υq

κp

Ψp ,

Υp

and it has the property that it is unaffected by the in-phase distortion (described in Section 11.1 below) which is recognised as common in MT data. Note, however, that the computation of phase tensor values may encounter difficulties for high condition numbers and singularities in **Zp** and **Zq**, which can be caused by strong distortion.

#### **8.1 Mohr diagram for the phase tensor**

The phase tensor, a 2 x 2 matrix, can also be represented by a Mohr diagram, as described by Weaver et al. (2003; 2006). The general form is shown in Fig. 9, following the convention for axes of Fig. 2a. Equations 6 to 15 can be adapted to apply to the quantities shown in Fig. 9. 14 Will-be-set-by-IN-TECH

κp,q

Fig. 7. The condition number *κp*,*<sup>q</sup>* displayed as a function of the anisotropy angle *λp*,*<sup>q</sup>* in

of the Mohr representations in Figs 2, 3 and 6, from Equation 41 the condition numbers are

*<sup>p</sup>*,*<sup>q</sup>* <sup>+</sup> *Cp*,*q*)/|(*Z<sup>L</sup>*

remembering Equation 16. Fig. 7 shows *κp*,*<sup>q</sup>* as a function of *λp*,*q*. While *λp*,*<sup>q</sup>* itself is a measure of condition, Fig. 7 demonstrates that *κp*,*<sup>q</sup>* is more sensitive than *λp*,*<sup>q</sup>* as ill-condition

An example is given in Fig. 8 of condition numbers determined for a set of data from a recent MT site, NQ142, in north Queensland, Australia. An increase in condition number above 10 monitors the decrease of the lesser principle value Ψ*p*,*q*, which becomes negligible. The variation of angle *θep*,*<sup>q</sup>* with period (T) stabilizes as condition numbers rise, an effect discussed

Caldwell et al. (2004) introduced a "phase tensor" based on the magnetotelluric tensor, and the concept was further developed by Bibby et al. (2005). The phase tensor is a real matrix **Φ**,

**Φ** = **Z**−**<sup>1</sup>**

and it has the property that it is unaffected by the in-phase distortion (described in Section 11.1 below) which is recognised as common in MT data. Note, however, that the computation of phase tensor values may encounter difficulties for high condition numbers and singularities

The phase tensor, a 2 x 2 matrix, can also be represented by a Mohr diagram, as described by Weaver et al. (2003; 2006). The general form is shown in Fig. 9, following the convention for axes of Fig. 2a. Equations 6 to 15 can be adapted to apply to the quantities shown in Fig. 9.

0 10 20 30 40 50 60 70 80 90 λp,q

*κp*,*<sup>q</sup>* = 2/(csc*λp*,*<sup>q</sup>* − 1) + 1 (43)

*<sup>p</sup>*,*<sup>q</sup>* − *Cp*,*q*)| (42)

**<sup>p</sup> Zq** (44)

is approached. Note that *κp*,*<sup>q</sup>* ≥ 1.

below in Section 11.

**8. The phase tensor**

defined by

given by

that is

0 10 20 30 40 50 60 70 80 90 λp,q

degrees, using linear (left) and logarithmic (right) scales for *κp*,*q*.

in **Zp** and **Zq**, which can be caused by strong distortion.

**8.1 Mohr diagram for the phase tensor**

*κp*,*<sup>q</sup>* = (*Z<sup>L</sup>*

κp,q

Fig. 8. The NQ142 data plotted as Mohr diagrams and analysed by SVD to give: angles *θep*, *θeq*, *θhp* and *θhq*; the principal values Υ*p*, Ψ*p*, Υ*<sup>q</sup>* and Ψ*q*; and the condition numbers *κ<sup>p</sup>* and *κq*. The variation of period with colour in the circles is the same as for the other plots. The Mohr diagrams follow the form of Fig. 2a, and have their impedance values multiplied by the square root of period (T) to make the plots more compact. Where circles enclose the origin, their lesser principal value (Ψ) is given an artificially high value to flag this circumstance; however condition numbers are computed nevertheless.

**8.2 Mohr diagram for the 2D phase tensor**

**8.3 Mohr diagram for the 1D phase tensor**

the phase value of the 1D MT response.

quadrant".

axes.

quadrant when both *Z*�

The condition for the phase of *Z*�

of their vertical axes, the *Z*�

**9. Invariants for the general phase tensor**

For the 2D case, the Mohr diagram reduces to a circle with centre on the Φ�

For the 1D case, the Mohr diagram reduces to a point on the Φ�

which neatly summarize 1D, 2D and 3D characteristics, respectively.

Adopting, in the rotated frame, the common definition for phase *φ*�

*xy p*

Phases greater than 90◦ are then possible for *Z*�

and *μp*,*<sup>q</sup>* must together be greater than 90◦, i.e.

and *Z*�

*φ*�

**10. Some complications illustrated by Mohr diagrams**

**10.1 Conditions for** *Zxy* **and** *Zyx* **phases out of quadrant**

given, are eigenvalues of the phase tensor. The two eigenvectors are orthogonal.

*TM* phase values are then the arctangents of the intercepts which, from the analysis already

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 97

To characterize the phase tensor just three invariants of rotation are needed, as pointed out by Caldwell et al. (2004), together with an angle defining the direction of the original observing axes. Options for invariants are evident from an inspection of Fig. 9, and include those adopted in Section 4 for Fig. 3. Weaver et al. (2003) chose *J*1, *J*<sup>2</sup> and *J*<sup>3</sup> as shown in Fig. 9,

For simple cases of induction in 1D layered media, *Zxy* phases will be in the first quadrant, and *Zyx* phases in the third quadrant. However, it is not uncommon with complicated geologic structure to find phases which are out of these expected quadrants, for some orientation of the measuring axes. Experimental or computational error may be invoked in explanation, when in fact a common cause is simply distortion. Such distortion may be understood by reference to Fig. 3. Clearly once, due say to distortion, angle *μ<sup>p</sup>* (or *μq*) increases to the point where the in-phase (or quadrature) circle first touches and then crosses the vertical *Z*�

for an appropriate direction of the observing axes the phase observed of *Zxy* will be "out of

quadrature parts should be negative. On the Mohr circle diagram, this condition is equivalent to either or both of the in-phase and quadrature circles for the data crossing the vertical axis.

With reference to Fig. 3, for the crossing of the vertical axis to occur, it can be seen that *λp*,*<sup>q</sup>*

*xyq* /*Z*� *xy p*

*xy* = arctan(*Z*�

where the signs of numerator and denominator are taken into account, then *φ*�

*xyq* are positive.

*xx* axis. The *TE* and

*xx* axis,

*xy* is in the first

*xx* axis, marking the tangent of

*xy* of

*xy*. If the circles however remain to the right

*xy* to go out of quadrant is that one or both of its in-phase and

*λp*,*<sup>q</sup>* + |*μp*,*q*| > 90◦ (46)

*xy* phase will be in quadrant for any orientation of the measuring

) (45)

Fig. 9. Diagram showing the Mohr representation of the general phase tensor. Axes for Φ� *yy* and Φ� *yx* could be added (following Fig. 2a). The principal values (equivalent to Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* in Fig. 6) are given by (*J*<sup>1</sup> <sup>2</sup> + *J*<sup>3</sup> <sup>2</sup>)1/2 <sup>±</sup> *<sup>J</sup>*2, and are denoted <sup>Φ</sup>*max* and <sup>Φ</sup>*min* by Caldwell et al. (2004).

The angle *γ*, which is a measure of 3D effects, is the seventh invariant of Weaver et al. (2000), described in Section 4.

The regional 2D strike direction chosen on phase considerations, as advocated by Bahr (1988) for example, is given by the point of maximum Φ� *xx*, which is the highest point of the circle as drawn. However when the regional structure is recognised as being 3D, the best estimate of geologic strike is recommended by Caldwell et al. (2004) as the orientation of the principal axes of the phase tensor, shown in Fig. 9 as "strike direction based on Φ*max*".

Note that a Mohr diagram for the phase tensor shows the variation of phase-tensor values as the observing axes are rotated, but that these values are generally not those of the usual MT *Zxy* phase.

Also note that a *Zxy* phase going "out of quadrant" does not imply that a principal value of **Φ** is negative. In fact det**Φ** is expected to be never negative, consistent with the principal values of **Φ** (Φ*max* and Φ*min*) being never negative.

Examples of Mohr diagrams for phase tensors, which illustrate some of these points, are given in Lilley & Weaver (2010).

#### **8.2 Mohr diagram for the 2D phase tensor**

16 Will-be-set-by-IN-TECH

Φ' xy

xx

*xx*, which is the highest point of the circle

*yy*

Φ' xx

(Φxy- Φyx)/2

J3

marks observed point

*yx* could be added (following Fig. 2a). The principal values (equivalent to Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* in

Fig. 9. Diagram showing the Mohr representation of the general phase tensor. Axes for Φ�

The angle *γ*, which is a measure of 3D effects, is the seventh invariant of Weaver et al. (2000),

The regional 2D strike direction chosen on phase considerations, as advocated by Bahr (1988)

as drawn. However when the regional structure is recognised as being 3D, the best estimate of geologic strike is recommended by Caldwell et al. (2004) as the orientation of the principal

Note that a Mohr diagram for the phase tensor shows the variation of phase-tensor values as the observing axes are rotated, but that these values are generally not those of the usual MT

Also note that a *Zxy* phase going "out of quadrant" does not imply that a principal value of **Φ** is negative. In fact det**Φ** is expected to be never negative, consistent with the principal values

Examples of Mohr diagrams for phase tensors, which illustrate some of these points, are given

axes of the phase tensor, shown in Fig. 9 as "strike direction based on Φ*max*".

marks general point after axes rotation

strike direction based on maximum Φ'

<sup>2</sup>)1/2 <sup>±</sup> *<sup>J</sup>*2, and are denoted <sup>Φ</sup>*max* and <sup>Φ</sup>*min* by Caldwell et al.

strike direction based on Φmax

γ

2θ'

J2

(Φxx+ Φyy)/2

0

<sup>2</sup> + *J*<sup>3</sup>

for example, is given by the point of maximum Φ�

of **Φ** (Φ*max* and Φ*min*) being never negative.

in Lilley & Weaver (2010).

and Φ�

(2004).

*Zxy* phase.

Fig. 6) are given by (*J*<sup>1</sup>

described in Section 4.

J1

For the 2D case, the Mohr diagram reduces to a circle with centre on the Φ� *xx* axis. The *TE* and *TM* phase values are then the arctangents of the intercepts which, from the analysis already given, are eigenvalues of the phase tensor. The two eigenvectors are orthogonal.

#### **8.3 Mohr diagram for the 1D phase tensor**

For the 1D case, the Mohr diagram reduces to a point on the Φ� *xx* axis, marking the tangent of the phase value of the 1D MT response.

#### **9. Invariants for the general phase tensor**

To characterize the phase tensor just three invariants of rotation are needed, as pointed out by Caldwell et al. (2004), together with an angle defining the direction of the original observing axes. Options for invariants are evident from an inspection of Fig. 9, and include those adopted in Section 4 for Fig. 3. Weaver et al. (2003) chose *J*1, *J*<sup>2</sup> and *J*<sup>3</sup> as shown in Fig. 9, which neatly summarize 1D, 2D and 3D characteristics, respectively.

## **10. Some complications illustrated by Mohr diagrams**

## **10.1 Conditions for** *Zxy* **and** *Zyx* **phases out of quadrant**

For simple cases of induction in 1D layered media, *Zxy* phases will be in the first quadrant, and *Zyx* phases in the third quadrant. However, it is not uncommon with complicated geologic structure to find phases which are out of these expected quadrants, for some orientation of the measuring axes. Experimental or computational error may be invoked in explanation, when in fact a common cause is simply distortion. Such distortion may be understood by reference to Fig. 3. Clearly once, due say to distortion, angle *μ<sup>p</sup>* (or *μq*) increases to the point where the in-phase (or quadrature) circle first touches and then crosses the vertical *Z*� *xx* axis, for an appropriate direction of the observing axes the phase observed of *Zxy* will be "out of quadrant".

Adopting, in the rotated frame, the common definition for phase *φ*� *xy* of

$$\phi\_{xy}^{\prime} = \arctan(Z\_{xy\_q}^{\prime}/Z\_{xy\_p}^{\prime}) \tag{45}$$

where the signs of numerator and denominator are taken into account, then *φ*� *xy* is in the first quadrant when both *Z*� *xy p* and *Z*� *xyq* are positive.

The condition for the phase of *Z*� *xy* to go out of quadrant is that one or both of its in-phase and quadrature parts should be negative. On the Mohr circle diagram, this condition is equivalent to either or both of the in-phase and quadrature circles for the data crossing the vertical axis. Phases greater than 90◦ are then possible for *Z*� *xy*. If the circles however remain to the right of their vertical axes, the *Z*� *xy* phase will be in quadrant for any orientation of the measuring axes.

With reference to Fig. 3, for the crossing of the vertical axis to occur, it can be seen that *λp*,*<sup>q</sup>* and *μp*,*<sup>q</sup>* must together be greater than 90◦, i.e.

$$
\lambda\_{p,q} + |\mu\_{p,q}| > 90^{\circ} \tag{46}
$$

and with distortion **d** present, affecting the measured electric field **E<sup>m</sup>** only,

**E<sup>m</sup>** = **ZmH** = **dZ<sup>r</sup>**

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 99

<sup>12</sup>; *<sup>Z</sup><sup>r</sup>*

<sup>12</sup>; *<sup>Z</sup><sup>r</sup>*

If the observation axes are north-south and east-west, say, and need to be rotated angle *θ* for

 0 *Z<sup>r</sup>* 12

> 0 *Z<sup>r</sup>* 12

*Zr* <sup>21</sup> 0

In co-ordinates aligned with the surficial geology the distortion amounts to scaling the electric fields by different amounts, *g*<sup>1</sup> and *g*2, in directions parallel and perpendicular to the surficial

> *g*<sup>1</sup> 0 0 *g*<sup>2</sup>

<sup>21</sup> sin *<sup>α</sup><sup>s</sup> <sup>g</sup>*1*Z<sup>r</sup>*

<sup>21</sup> cos *<sup>α</sup><sup>s</sup>* <sup>−</sup>*g*2*Z<sup>r</sup>*

At this stage the in-phase and quadrature parts can be considered separately. Avoiding the encumbrance of additional notation, now assume that just the in-phase part of the measured *E* field is being considered, giving information on just the in-phase parts of the regional

Then following the form of Equation 32, SVD can be carried out on the central tensor in

= **R**(−*ηe*)

<sup>12</sup> cos *α<sup>s</sup>*

<sup>12</sup> sin *α<sup>s</sup>*

 0 Υ −Ψ 0

= **R**(−*αs*)

where due to the restricted distortion model the *d*<sup>21</sup> element has the value *d*12, *g*<sup>1</sup> and *g*<sup>2</sup> are real constants (frequency independent), and *α<sup>s</sup>* is the angle from the regional strike coordinates to the strike coordinates of the surficial structure causing the distortion. The

 **R**(*θ*) *Hx Hy* 

 **R**(*θ*) *Hx Hy* 

*Zr* <sup>21</sup> 0

For the case of 2D regional structure, with measuring axes aligned with regional strike,

**Z˜ <sup>r</sup>** = [0, *Z<sup>r</sup>*

**E˜ <sup>m</sup>** = **d**[0, *Z<sup>r</sup>*

where the overscore**˜**indicates observing axes aligned with regional strike.

alignment with regional strike, then Equation 51 (without distortion) gives

**R**(*θ*) *E<sup>m</sup> x E<sup>m</sup> y* =

**R**(*θ*) *E<sup>m</sup> x E<sup>m</sup> y* = **d**

> *d*<sup>11</sup> *d*<sup>12</sup> *<sup>d</sup>*<sup>12</sup> *<sup>d</sup>*<sup>22</sup>

decomposition of **d** in Equation 54 is displayed in Fig. 10.

**<sup>E</sup><sup>m</sup>** <sup>=</sup> **<sup>R</sup>**(−*<sup>θ</sup>* <sup>−</sup> *<sup>α</sup>s*)

<sup>21</sup> sin *<sup>α</sup><sup>s</sup> <sup>g</sup>*1*Z<sup>r</sup>*

<sup>21</sup> cos *<sup>α</sup><sup>s</sup>* <sup>−</sup>*g*2*Z<sup>r</sup>*

Substituting the expression for **d** from Equation 54 into Equation 53 gives

*g*1*Z<sup>r</sup>*

*g*2*Z<sup>r</sup>*

<sup>12</sup> cos *α<sup>s</sup>*

<sup>12</sup> sin *α<sup>s</sup>*

strike (Zhang et al., 1987). Distortion matrix **d** then has the form

and

and so

and with distortion

impedance tensor **Zr**.

Equation 55, expanding it as

*g*1*Z<sup>r</sup>*

*g*2*Z<sup>r</sup>*

**E<sup>m</sup>** = **dE<sup>r</sup>** (48)

**H** (49)

<sup>21</sup>, 0] (50)

<sup>21</sup>, 0]**H** (51)

**R**(*αs*) (54)

**R**(*θ*)**H** (55)

**R**(*ηh*) (56)

(52)

(53)

With reference to the *Z*� *yy p* and *Z*� *yx <sup>p</sup>* axes also included in Fig. 2a, consequences for the phases of other elements are also clear. The symmetry of the figure is such that if the circle crosses the *Z*� *xx <sup>p</sup>* axis to the left, it will also cross the *Z*� *yy <sup>p</sup>* axis to the right. Thus if phases out of quadrant are possible for *Z*� *xy*, so are they also possible for *Z*� *yx*. However, the latter will occur for rotations 90◦ different from the former.

Similar considerations will give rules regarding the quadrants to be expected for *Z*� *xx* and *Z*� *yy* phases. An analysis of Fig. 2a shows that these phases will generally change quadrant with a rotation of the observing axes. An exception, when they would not do so, would be if the circle in Fig. 2a did not cross the horizontal *Z*� *xy <sup>p</sup>* axis, but stayed completely above (or below) this axis (and the circle for the quadrature part of the tensor behaved similarly).

#### **10.2 Strike coordinates from phase considerations**

The situation of *Zxy* phase out of quadrant (for example for the part of the circle to the left of the vertical axis in Fig. 5c) causes difficulty in the symmetric-antisymmetric and the symmetric-2D partitions of Sections 5.1 and 5.2. First, if the circle in Fig. 5c is simply moved down so that its centre lies on the horizontal axis, the left-hand intercept with the horizontal axis is less than zero. The circle has enclosed the origin, and a negative value is obtained for the lesser of *TE* and *TM*, contrary to expectation (see Section 13.1).

Similarly there are methods which seek, as geologic strike coordinates, those given by maximum and minimum *Z*� *xy* phase values when the axes are rotated. These methods will find erratic phase behaviour when either or both of the in-phase and quadrature circles cross the vertical axis as in Fig. 5c.

However *Z*� *xy* phase will usually be at or near a maximum (or minimum) at the right-hand side of a circle which on its left-hand side crosses the vertical axis. Thus sensible results for maximum or minimum phase can be expected, if based only on a determination from the right-hand side.

## **11. Geologic interpretation of angles found by SVD, when condition numbers are high**

Singular value decomposition of observed magnetotelluric tensors, taking the in-phase and quadrature parts separately, often produces directions which are frequency independent below a certain frequency which is found to be common to both in-phase and quadrature parts. This section examines how such directions may be related to the angles which arise in basic distortion models of MT data. It is found that while MT data are expected in general to be frequency dependent, there are particular limiting cases which are constant with frequency. Understanding such observed data may be useful in interpretation.

#### **11.1 The Smith (1995) model for local distortion**

This section applies the Smith (1995) description of local "static" distortion, which follows Bahr (1988) and Groom & Bailey (1989). In the absence of local distortion the measured magnetic field **H** is related to the regional electric field **E<sup>r</sup>** by the tensor **Zr**:

$$\mathbf{E}^{\mathbf{r}} = \mathbf{Z}^{\mathbf{r}} \mathbf{H} \tag{47}$$

and with distortion **d** present, affecting the measured electric field **E<sup>m</sup>** only,

$$\mathbf{E}^{\mathbf{m}} = \mathbf{d} \mathbf{E}^{\mathbf{r}} \tag{48}$$

and

18 Will-be-set-by-IN-TECH

of other elements are also clear. The symmetry of the figure is such that if the circle crosses

phases. An analysis of Fig. 2a shows that these phases will generally change quadrant with a rotation of the observing axes. An exception, when they would not do so, would be if the

The situation of *Zxy* phase out of quadrant (for example for the part of the circle to the left of the vertical axis in Fig. 5c) causes difficulty in the symmetric-antisymmetric and the symmetric-2D partitions of Sections 5.1 and 5.2. First, if the circle in Fig. 5c is simply moved down so that its centre lies on the horizontal axis, the left-hand intercept with the horizontal axis is less than zero. The circle has enclosed the origin, and a negative value is obtained for

Similarly there are methods which seek, as geologic strike coordinates, those given by

find erratic phase behaviour when either or both of the in-phase and quadrature circles cross

side of a circle which on its left-hand side crosses the vertical axis. Thus sensible results for maximum or minimum phase can be expected, if based only on a determination from the

**11. Geologic interpretation of angles found by SVD, when condition numbers are**

Singular value decomposition of observed magnetotelluric tensors, taking the in-phase and quadrature parts separately, often produces directions which are frequency independent below a certain frequency which is found to be common to both in-phase and quadrature parts. This section examines how such directions may be related to the angles which arise in basic distortion models of MT data. It is found that while MT data are expected in general to be frequency dependent, there are particular limiting cases which are constant with frequency.

This section applies the Smith (1995) description of local "static" distortion, which follows Bahr (1988) and Groom & Bailey (1989). In the absence of local distortion the measured

**E<sup>r</sup>** = **Z<sup>r</sup>**

*xy* phase will usually be at or near a maximum (or minimum) at the right-hand

*xy*, so are they also possible for *Z*�

Similar considerations will give rules regarding the quadrants to be expected for *Z*�

this axis (and the circle for the quadrature part of the tensor behaved similarly).

*yx <sup>p</sup>* axes also included in Fig. 2a, consequences for the phases

*xy* phase values when the axes are rotated. These methods will

*yy <sup>p</sup>* axis to the right. Thus if phases out of

*xy <sup>p</sup>* axis, but stayed completely above (or below)

**H** (47)

*yx*. However, the latter will occur

*xx* and *Z*�

*yy*

With reference to the *Z*�

quadrant are possible for *Z*�

maximum and minimum *Z*�

the vertical axis as in Fig. 5c.

However *Z*�

right-hand side.

**high**

the *Z*�

*yy p*

for rotations 90◦ different from the former.

circle in Fig. 2a did not cross the horizontal *Z*�

**10.2 Strike coordinates from phase considerations**

the lesser of *TE* and *TM*, contrary to expectation (see Section 13.1).

Understanding such observed data may be useful in interpretation.

magnetic field **H** is related to the regional electric field **E<sup>r</sup>** by the tensor **Zr**:

**11.1 The Smith (1995) model for local distortion**

and *Z*�

*xx <sup>p</sup>* axis to the left, it will also cross the *Z*�

$$\mathbf{E}^{\mathbf{m}} = \mathbf{Z}^{\mathbf{m}} \mathbf{H} = \mathbf{d} \mathbf{Z}^{\mathbf{r}} \mathbf{H} \tag{49}$$

For the case of 2D regional structure, with measuring axes aligned with regional strike,

$$\mathbf{Z}^{\mathbf{r}} = [0, Z\_{12}^{r}; Z\_{21}^{r}, 0] \tag{50}$$

and so

$$\mathbf{\tilde{E}^m} = \mathbf{d}[0, Z\_{12}^r; Z\_{21}^r, 0] \mathbf{H} \tag{51}$$

where the overscore**˜**indicates observing axes aligned with regional strike.

If the observation axes are north-south and east-west, say, and need to be rotated angle *θ* for alignment with regional strike, then Equation 51 (without distortion) gives

$$\mathbf{R}(\theta) \begin{bmatrix} E\_x^m \\ E\_y^m \end{bmatrix} = \begin{bmatrix} 0 & Z\_{12}^r \\ Z\_{21}^r & 0 \end{bmatrix} \mathbf{R}(\theta) \begin{bmatrix} H\_x \\ H\_y \end{bmatrix} \tag{52}$$

and with distortion

$$\mathbf{R}(\theta) \begin{bmatrix} E\_x^m \\ E\_y^m \end{bmatrix} = \mathbf{d} \begin{bmatrix} 0 & Z\_{12}^r \\ Z\_{21}^r & 0 \end{bmatrix} \mathbf{R}(\theta) \begin{bmatrix} H\_x \\ H\_y \end{bmatrix} \tag{53}$$

In co-ordinates aligned with the surficial geology the distortion amounts to scaling the electric fields by different amounts, *g*<sup>1</sup> and *g*2, in directions parallel and perpendicular to the surficial strike (Zhang et al., 1987). Distortion matrix **d** then has the form

$$
\begin{bmatrix} d\_{11} \ d\_{12} \\ d\_{12} \ d\_{22} \end{bmatrix} = \mathbf{R}(-a\_s) \begin{bmatrix} \mathcal{g}\_1 & 0 \\ 0 & \mathcal{g}\_2 \end{bmatrix} \mathbf{R}(a\_s) \tag{54}
$$

where due to the restricted distortion model the *d*<sup>21</sup> element has the value *d*12, *g*<sup>1</sup> and *g*<sup>2</sup> are real constants (frequency independent), and *α<sup>s</sup>* is the angle from the regional strike coordinates to the strike coordinates of the surficial structure causing the distortion. The decomposition of **d** in Equation 54 is displayed in Fig. 10.

Substituting the expression for **d** from Equation 54 into Equation 53 gives

$$\mathbf{E}^{\mathbf{m}} = \mathbf{R}(-\theta - \mathfrak{a}\_{\mathrm{s}}) \begin{bmatrix} g\_1 Z\_{21}^{\mathrm{r}} \sin \mathfrak{a}\_{\mathrm{s}} & g\_1 Z\_{12}^{\mathrm{r}} \cos \mathfrak{a}\_{\mathrm{s}}\\ g\_2 Z\_{21}^{\mathrm{r}} \cos \mathfrak{a}\_{\mathrm{s}} & -g\_2 Z\_{12}^{\mathrm{r}} \sin \mathfrak{a}\_{\mathrm{s}} \end{bmatrix} \mathbf{R}(\theta) \mathbf{H} \tag{55}$$

At this stage the in-phase and quadrature parts can be considered separately. Avoiding the encumbrance of additional notation, now assume that just the in-phase part of the measured *E* field is being considered, giving information on just the in-phase parts of the regional impedance tensor **Zr**.

Then following the form of Equation 32, SVD can be carried out on the central tensor in Equation 55, expanding it as

$$
\begin{bmatrix}
\text{g}\_1 Z\_{21}^{\prime} \sin \mathfrak{a}\_s & \text{g}\_1 Z\_{12}^{\prime} \cos \mathfrak{a}\_s \\
\text{g}\_2 Z\_{21}^{\prime} \cos \mathfrak{a}\_s & -\text{g}\_2 Z\_{12}^{\prime} \sin \mathfrak{a}\_s
\end{bmatrix} = \mathbf{R}(-\eta\_\ell) \begin{bmatrix}
\text{0} & \text{Y} \\
\end{bmatrix} \mathbf{R}(\eta\_h) \tag{56}
$$

**11.2 First limiting case**

and

and

strike.

case.

and

and

strike.

Further, note that if *Z<sup>r</sup>*

**11.3 Second limiting case**

Now consider, in Equations 57, 58, 59 and 60, that *g*<sup>2</sup> � *g*1, so that terms in *g*<sup>1</sup> are ignored

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 101

<sup>12</sup> sin *<sup>α</sup>s*)<sup>2</sup> + (*Z<sup>r</sup>*

*<sup>η</sup><sup>h</sup>* <sup>=</sup> <sup>−</sup> arctan(*Z<sup>r</sup>*

<sup>12</sup>, *<sup>Z</sup><sup>r</sup>*

simple result of zero. Fed back into Equation 61 it is seen that SVD of an MT tensor now gives, in the equivalent of *θe* in Equation 31, a result for (*θ* + *αs*): that is, the direction of surficial 2D

with (or normal to) the surficial strike. This case can be seen to be that where the regional 2D anisotropy is out-weighed by the surficial anisotropy, so that the former approximates a 1D

A second limiting case might be gross inequality in the *TE* and *TM* components of the regional

be ignored with respect to terms involving the former. Then Equations 57, 58, 59 and 60 give

Fed back into Equation 61, it is seen that SVD of an MT tensor for the case where *η<sup>h</sup>* = 0 now gives, in the equivalent of *θh* in Equation 31, a result for *θ*: that is, the direction of regional 2D

The two limiting cases, where realised, give angles of interest. Also, it should be noted that

12|�|*Z<sup>r</sup>*

<sup>12</sup>[(*g*<sup>2</sup> sin *<sup>α</sup>s*)<sup>2</sup> + (*g*<sup>1</sup> cos *<sup>α</sup>s*)2]

<sup>21</sup> cos *<sup>α</sup>s*)2]

<sup>12</sup>tan *<sup>α</sup>s*/*Z<sup>r</sup>*

<sup>21</sup> = −1, then *η<sup>h</sup>* = *αs*, and the magnetic axes are also aligned

1/2 (62)

<sup>21</sup>) (65)

Ψ = 0 (63)

*η<sup>e</sup>* = 0 (64)

<sup>21</sup> and *α<sup>s</sup>* still occur together, for *η<sup>e</sup>* there is the

<sup>21</sup>|, so that terms involving the latter may

Ψ = 0 (67)

*η<sup>h</sup>* = 0 (69)

*η<sup>e</sup>* = − arctan(*g*2tan *αs*/*g*1) (68)

1/2 (66)

with respect to terms in *g*2. The results are then obtained that

For the angles *η<sup>e</sup>* and *η<sup>h</sup>* the limiting case gives

While in Equation 65 for *η<sup>h</sup>* the values *Z<sup>r</sup>*

12/*Z<sup>r</sup>*

2D impedance. Consider the case where <sup>|</sup>*Z<sup>r</sup>*

For the angles *η<sup>e</sup>* and *η<sup>h</sup>* this limiting case gives

**12. Discussion and example**

Υ = *Z<sup>r</sup>*

both cases correspond to the MT tensor being singular, or nearly so.

Υ = *g*2[(*Z<sup>r</sup>*

Fig. 10. Mohr diagrams for the decomposition in Equation 54. On the left the diagram shows the distortion matrix **d**. It is positive definite, due to the assumption of strictly 2D surficial anisotropy, and so has orthogonal eigenvectors. The radius of the circle is (*g*<sup>2</sup> − *g*1)/2. The numerical matrix taken for the example is [3.5, 2.6; 2.6, 6.5], and the directions and magnitudes of its eigenvectors are shown by rotating the radial arm until it is parallel to the vertical axis. The eigenvalues then are *g*<sup>1</sup> = 2 and *g*<sup>2</sup> = 8. The angle 2*α<sup>s</sup>* on the figure is 60◦, indicating that a rotation of axes anticlockwise through 30◦ is necessary to change the matrix to [2, 0: 0, 8]. On the right is shown the matrix decomposition as in Equation 54. The three diagrams represent the three matrices on the right-hand side of that equation, in turn. The first Mohr diagram represents a rotation by angle (−*αs*). The second Mohr diagram is symmetric, representing a positive definite matrix according to the scaling factors *g*<sup>1</sup> and *g*2. The third matrix is like the first, representing a rotation by angle (+*αs*). Note the changes of scale between the diagrams.

where the angles *η<sup>e</sup>* and *ηh*, and the principal values Υ and Ψ are given by

$$\eta\_{\ell} + \eta\_{h} = \arctan\left[\frac{-g\_2 Z\_{12}^r \sin \alpha\_s - g\_1 Z\_{21}^r \sin \alpha\_s}{g\_1 Z\_{12}^r \cos \alpha\_s + g\_2 Z\_{21}^r \cos \alpha\_s}\right] \tag{57}$$

$$\eta\_{\ell} - \eta\_{h} = \arctan\left[\frac{-g\_{2}Z\_{12}^{r}\sin\alpha\_{s} + g\_{1}Z\_{21}^{r}\sin\alpha\_{s}}{g\_{1}Z\_{12}^{r}\cos\alpha\_{s} - g\_{2}Z\_{21}^{r}\cos\alpha\_{s}}\right] \tag{58}$$

$$\mathbf{Y} + \mathbf{Y} = \left[ (\mathbf{g}\_1 \mathbf{Z}\_{21}^r \sin \mathfrak{a}\_s - \mathbf{g}\_2 \mathbf{Z}\_{12}^r \sin \mathfrak{a}\_s)^2 + (\mathbf{g}\_2 \mathbf{Z}\_{21}^r \cos \mathfrak{a}\_s - \mathbf{g}\_1 \mathbf{Z}\_{12}^r \cos \mathfrak{a}\_s)^2 \right]^{1/2} \tag{59}$$

and

$$\Psi - \Psi = \left[ (\mathcal{g}\_1 \mathcal{Z}\_{21}^r \sin \mathfrak{a}\_s + \mathcal{g}\_2 \mathcal{Z}\_{12}^r \sin \mathfrak{a}\_s)^2 + (\mathcal{g}\_1 \mathcal{Z}\_{12}^r \cos \mathfrak{a}\_s + \mathcal{g}\_2 \mathcal{Z}\_{21}^r \cos \mathfrak{a}\_s)^2 \right]^{1/2} \tag{60}$$

Replacing into Equation 55 the matrix in its expanded form as in Equation 56, and combining adjoining rotation matrices, gives

$$\mathbf{E}^{\mathbf{m}} = \mathbf{R}(-\theta - \mathfrak{a}\_s - \eta\_\varepsilon) \begin{bmatrix} 0 & \mathbf{Y} \\ -\mathbf{y} & 0 \end{bmatrix} \mathbf{R}(\theta + \eta\_h) \mathbf{H} \tag{61}$$

Thus in the SVD of an MT tensor resulting from 2D surficial distortion of a regional 2D structure, the regional strike *θ* and the various distortion quantities *g*1, *g*<sup>2</sup> and *α<sup>s</sup>* occur in a way which does not allow their straightforward individual solution, because the SVD produces values for (*θ* + *α<sup>s</sup>* + *ηe*) and (*θ* + *ηh*). It is of interest however to examine two limiting cases, as in the following sections.

#### **11.2 First limiting case**

Now consider, in Equations 57, 58, 59 and 60, that *g*<sup>2</sup> � *g*1, so that terms in *g*<sup>1</sup> are ignored with respect to terms in *g*2. The results are then obtained that

$$\mathbf{Y} = \mathbf{g}\_2[(Z\_{12}^r \sin \mathfrak{a}\_s)^2 + (Z\_{21}^r \cos \mathfrak{a}\_s)^2]^{1/2} \tag{62}$$

and

20 Will-be-set-by-IN-TECH

cos α<sup>s</sup> cos α<sup>s</sup> <sup>2</sup>α<sup>s</sup>

0 0 0 0

Fig. 10. Mohr diagrams for the decomposition in Equation 54. On the left the diagram shows the distortion matrix **d**. It is positive definite, due to the assumption of strictly 2D surficial anisotropy, and so has orthogonal eigenvectors. The radius of the circle is (*g*<sup>2</sup> − *g*1)/2. The

magnitudes of its eigenvectors are shown by rotating the radial arm until it is parallel to the vertical axis. The eigenvalues then are *g*<sup>1</sup> = 2 and *g*<sup>2</sup> = 8. The angle 2*α<sup>s</sup>* on the figure is 60◦, indicating that a rotation of axes anticlockwise through 30◦ is necessary to change the matrix to [2, 0: 0, 8]. On the right is shown the matrix decomposition as in Equation 54. The three diagrams represent the three matrices on the right-hand side of that equation, in turn. The first Mohr diagram represents a rotation by angle (−*αs*). The second Mohr diagram is symmetric, representing a positive definite matrix according to the scaling factors *g*<sup>1</sup> and *g*2. The third matrix is like the first, representing a rotation by angle (+*αs*). Note the changes of


2 2

d'11

d'12

scale between the diagrams.

Υ + Ψ = [(*g*1*Z<sup>r</sup>*

<sup>Υ</sup> <sup>−</sup> <sup>Ψ</sup> = [(*g*1*Z<sup>r</sup>*

adjoining rotation matrices, gives

cases, as in the following sections.

and

g1 g1

<sup>2</sup>**=** <sup>2</sup>

numerical matrix taken for the example is [3.5, 2.6; 2.6, 6.5], and the directions and

where the angles *η<sup>e</sup>* and *ηh*, and the principal values Υ and Ψ are given by

<sup>−</sup>*g*2*Z<sup>r</sup>*

*g*1*Z<sup>r</sup>*

<sup>−</sup>*g*2*Z<sup>r</sup>*

*g*1*Z<sup>r</sup>*

<sup>12</sup> sin *<sup>α</sup>s*)<sup>2</sup> + (*g*2*Z<sup>r</sup>*

<sup>12</sup> sin *<sup>α</sup>s*)<sup>2</sup> + (*g*1*Z<sup>r</sup>*

 0 Υ −Ψ 0

Replacing into Equation 55 the matrix in its expanded form as in Equation 56, and combining

Thus in the SVD of an MT tensor resulting from 2D surficial distortion of a regional 2D structure, the regional strike *θ* and the various distortion quantities *g*1, *g*<sup>2</sup> and *α<sup>s</sup>* occur in a way which does not allow their straightforward individual solution, because the SVD produces values for (*θ* + *α<sup>s</sup>* + *ηe*) and (*θ* + *ηh*). It is of interest however to examine two limiting

<sup>12</sup> sin *<sup>α</sup><sup>s</sup>* <sup>−</sup> *<sup>g</sup>*1*Z<sup>r</sup>*

<sup>12</sup> sin *<sup>α</sup><sup>s</sup>* <sup>+</sup> *<sup>g</sup>*1*Z<sup>r</sup>*

<sup>12</sup> cos *<sup>α</sup><sup>s</sup>* <sup>+</sup> *<sup>g</sup>*2*Z<sup>r</sup>*

<sup>12</sup> cos *<sup>α</sup><sup>s</sup>* <sup>−</sup> *<sup>g</sup>*2*Z<sup>r</sup>*

<sup>21</sup> sin *α<sup>s</sup>*

<sup>21</sup> sin *α<sup>s</sup>*

<sup>21</sup> cos *α<sup>s</sup>*

<sup>21</sup> cos *<sup>α</sup><sup>s</sup>* <sup>−</sup> *<sup>g</sup>*1*Z<sup>r</sup>*

<sup>12</sup> cos *<sup>α</sup><sup>s</sup>* <sup>+</sup> *<sup>g</sup>*2*Z<sup>r</sup>*

<sup>12</sup> cos *<sup>α</sup>s*)2]

<sup>21</sup> cos *<sup>α</sup>s*)2]

**R**(*θ* + *ηh*)**H** (61)

(57)

(58)

1/2 (59)

1/2 (60)

<sup>21</sup> cos *α<sup>s</sup>*

*η<sup>e</sup>* + *η<sup>h</sup>* = arctan

*η<sup>e</sup>* − *η<sup>h</sup>* = arctan

<sup>21</sup> sin *<sup>α</sup><sup>s</sup>* <sup>−</sup> *<sup>g</sup>*2*Z<sup>r</sup>*

<sup>21</sup> sin *<sup>α</sup><sup>s</sup>* <sup>+</sup> *<sup>g</sup>*2*Z<sup>r</sup>*

**<sup>E</sup><sup>m</sup>** <sup>=</sup> **<sup>R</sup>**(−*<sup>θ</sup>* <sup>−</sup> *<sup>α</sup><sup>s</sup>* <sup>−</sup> *<sup>η</sup>e*)

g2 g2

$$\mathbf{A} = \mathbf{0} \tag{63}$$

For the angles *η<sup>e</sup>* and *η<sup>h</sup>* the limiting case gives

$$
\eta\_{\ell} = 0 \tag{64}
$$

and

$$\eta\_h = -\arctan(Z\_{12}^r \tan \alpha\_s / Z\_{21}^r) \tag{65}$$

While in Equation 65 for *η<sup>h</sup>* the values *Z<sup>r</sup>* <sup>12</sup>, *<sup>Z</sup><sup>r</sup>* <sup>21</sup> and *α<sup>s</sup>* still occur together, for *η<sup>e</sup>* there is the simple result of zero. Fed back into Equation 61 it is seen that SVD of an MT tensor now gives, in the equivalent of *θe* in Equation 31, a result for (*θ* + *αs*): that is, the direction of surficial 2D strike.

Further, note that if *Z<sup>r</sup>* 12/*Z<sup>r</sup>* <sup>21</sup> = −1, then *η<sup>h</sup>* = *αs*, and the magnetic axes are also aligned with (or normal to) the surficial strike. This case can be seen to be that where the regional 2D anisotropy is out-weighed by the surficial anisotropy, so that the former approximates a 1D case.

#### **11.3 Second limiting case**

A second limiting case might be gross inequality in the *TE* and *TM* components of the regional 2D impedance. Consider the case where <sup>|</sup>*Z<sup>r</sup>* 12|�|*Z<sup>r</sup>* <sup>21</sup>|, so that terms involving the latter may be ignored with respect to terms involving the former. Then Equations 57, 58, 59 and 60 give

$$\mathbf{Y} = Z\_{12}^{r} [(\mathbf{g}\_2 \sin \alpha\_s)^2 + (\mathbf{g}\_1 \cos \alpha\_s)^2]^{1/2} \tag{66}$$

and

$$\mathbf{\hat{A}} = \mathbf{0} \tag{67}$$

For the angles *η<sup>e</sup>* and *η<sup>h</sup>* this limiting case gives

$$\eta\_{\ell} = -\arctan(g\_2 \tan \alpha\_s / g\_1) \tag{68}$$

and

$$
\eta\_h = 0 \tag{69}
$$

Fed back into Equation 61, it is seen that SVD of an MT tensor for the case where *η<sup>h</sup>* = 0 now gives, in the equivalent of *θh* in Equation 31, a result for *θ*: that is, the direction of regional 2D strike.

#### **12. Discussion and example**

The two limiting cases, where realised, give angles of interest. Also, it should be noted that both cases correspond to the MT tensor being singular, or nearly so.

observed point (Zxxp=-3, Zxyp =3)

Hx Hx'

θhp

v2

Hy

*<sup>p</sup>***Z***<sup>p</sup>* will each have a negative eigenvalue. However

**Z***<sup>p</sup>* = [−3, 3; −1, 5] (70)

Hy'

Z'xxp

1

Ex

θep

**13.3 Circles capturing the origin: a note on negative determinants**

*<sup>p</sup>* and **Z<sup>T</sup>**

To demonstrate this point, consider the hypothetical matrix

Ψ <sup>p</sup>

Ey'

origin defines positive).

negative, then the matrices **Z***p***Z<sup>T</sup>**

in Fig. 12.

0 4

2θep

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 103

θep - θhp

u1

<sup>v</sup> <sup>u</sup> <sup>1</sup> <sup>2</sup>

Ey

Fig. 12. The Mohr diagram for the matrix in Equation 70. The circle encloses the origin of axes because the matrix has a negative determinant. The values of *θep*, *θhp*, Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* are evident as 51◦, 25◦, 6.3 and −1.8 (Ψ*<sup>p</sup>* being of negative sign when the direction of Υ*<sup>p</sup>* from the

It is also observed that circles do not capture or enclose their origins of axes, except for reasons of obvious error, or usual errors associated with singularity. The author suggests that this observed behaviour corresponds to the prohibition of components of negative resistivity in the electrical conductivity structure. Again a proof regarding this possibility would be welcome. Negative determinants (no matter how they arise) need care in their analysis, as the following example illustrates. If, in contrast to the example in Section 6, the determinant (say of **Z***p*) is

in a conventional SVD these negative eigenvalues are taken as positive singular values, and to allow for this action the sign is also changed of the eigenvector of **U** which they multiply.

which has a negative determinant. The Mohr diagram representation of the matrix is shown

Ex'

Υ p

Z'xyp

marks "observed" points

Fig. 11. Hypothetical example of a 2D case with the in-phase and quadrature radial arms anti-parallel. Such cases appear to be prohibited in nature, but inadvertently can be posed as numerical examples.

The example in Fig. 8, for frequencies less that 10 Hz, gives *θep*,*<sup>q</sup>* values which are fixed at −67◦, while *θhp*,*<sup>q</sup>* values continue to vary with period. The model of Section 11.1 is therefore a good candidate for the interpretation of the data. Applying the results of Section 11.2, surficial strike is determined at −67◦ (±90◦ allowing for the usual ambiguity).

However the observation that there is a critical frequency, say *fc*, only below which *θep*,*<sup>q</sup>* values are frequency independent (and in the case of Fig. 8 are stable at −67◦), implies a contradiction to the initial distortion model in Section 11.1, which was frequency independent, generally. The existence of such a "critical frequency" suggests distortion which is "at a distance", rather than immediately local.

#### **13. Some remaining problems**

#### **13.1** *TE* **and** *TM* **modes both positive?**

It is well known that for the 2D case the *TM* mode is always positive (Weidelt & Kaikkonen, 1994), but proof has not been achieved for the *TE* case. It can be seen from Fig. 2b that were the *TE* case to be negative in either in-phase or quadrature part (or both) then the appropriate circles would enclose their origin of axes. Thus a general proof that circles cannot enclose their origins (discussed in Section 13.3 below) would also prove that the *TE* mode can never be negative.

#### **13.2 Radial arms anti-parallel**

The 2D example given in Fig. 2b, where the in-phase and quadrature radial arms are parallel, represents the common case observed without exception in the experience of the present author. However, in principle the case shown in Fig. 11 where the in-phase and quadrature radial arms are anti-parallel would also be two dimensional. If there is a theoretical reason why cases with radial arms anti-parallel as in Fig. 11 do not occur in practice it would be an advance to have it clarified.

22 Will-be-set-by-IN-TECH

Z'xxq

β

In-phase Quadrature

Z'xyp

β

strike is determined at −67◦ (±90◦ allowing for the usual ambiguity).

marks "observed" points

Fig. 11. Hypothetical example of a 2D case with the in-phase and quadrature radial arms anti-parallel. Such cases appear to be prohibited in nature, but inadvertently can be posed as

The example in Fig. 8, for frequencies less that 10 Hz, gives *θep*,*<sup>q</sup>* values which are fixed at −67◦, while *θhp*,*<sup>q</sup>* values continue to vary with period. The model of Section 11.1 is therefore a good candidate for the interpretation of the data. Applying the results of Section 11.2, surficial

However the observation that there is a critical frequency, say *fc*, only below which *θep*,*<sup>q</sup>* values are frequency independent (and in the case of Fig. 8 are stable at −67◦), implies a contradiction to the initial distortion model in Section 11.1, which was frequency independent, generally. The existence of such a "critical frequency" suggests distortion which is "at a

It is well known that for the 2D case the *TM* mode is always positive (Weidelt & Kaikkonen, 1994), but proof has not been achieved for the *TE* case. It can be seen from Fig. 2b that were the *TE* case to be negative in either in-phase or quadrature part (or both) then the appropriate circles would enclose their origin of axes. Thus a general proof that circles cannot enclose their origins (discussed in Section 13.3 below) would also prove that the *TE* mode can never

The 2D example given in Fig. 2b, where the in-phase and quadrature radial arms are parallel, represents the common case observed without exception in the experience of the present author. However, in principle the case shown in Fig. 11 where the in-phase and quadrature radial arms are anti-parallel would also be two dimensional. If there is a theoretical reason why cases with radial arms anti-parallel as in Fig. 11 do not occur in practice it would be an

Z'xyq <sup>0</sup> <sup>0</sup>

Z'xxp

numerical examples.

distance", rather than immediately local.

**13.1** *TE* **and** *TM* **modes both positive?**

**13. Some remaining problems**

**13.2 Radial arms anti-parallel**

advance to have it clarified.

be negative.

Fig. 12. The Mohr diagram for the matrix in Equation 70. The circle encloses the origin of axes because the matrix has a negative determinant. The values of *θep*, *θhp*, Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* are evident as 51◦, 25◦, 6.3 and −1.8 (Ψ*<sup>p</sup>* being of negative sign when the direction of Υ*<sup>p</sup>* from the origin defines positive).

#### **13.3 Circles capturing the origin: a note on negative determinants**

It is also observed that circles do not capture or enclose their origins of axes, except for reasons of obvious error, or usual errors associated with singularity. The author suggests that this observed behaviour corresponds to the prohibition of components of negative resistivity in the electrical conductivity structure. Again a proof regarding this possibility would be welcome.

Negative determinants (no matter how they arise) need care in their analysis, as the following example illustrates. If, in contrast to the example in Section 6, the determinant (say of **Z***p*) is negative, then the matrices **Z***p***Z<sup>T</sup>** *<sup>p</sup>* and **Z<sup>T</sup>** *<sup>p</sup>***Z***<sup>p</sup>* will each have a negative eigenvalue. However in a conventional SVD these negative eigenvalues are taken as positive singular values, and to allow for this action the sign is also changed of the eigenvector of **U** which they multiply.

To demonstrate this point, consider the hypothetical matrix

$$\mathbf{Z}\_p = [-\mathbf{3}, \mathbf{3}; -1, \mathbf{5}] \tag{70}$$

which has a negative determinant. The Mohr diagram representation of the matrix is shown in Fig. 12.

Invariants of an observed tensor under axes rotation are evident from inspection of such Mohr diagrams. The invariants may be used as indicators of the geologic dimensionality of the tensor; they may also be useful in the inversion and interpretation of MT data. Mohr diagrams of the MT phase tensor also show invariants which measure geologic dimensionality, and

Magnetotelluric Tensor Decomposition: Insights from Linear Algebra and Mohr Diagrams 105

Attention has been given to monitoring the singularity of an MT tensor with a condition number, as tensors which approach singularity are common. Carrying out SVD on observed data with high condition numbers can indicate a direction of surficial 2D strike, for the simple case of the surficial distortion of a regional 2D MT response. In some circumstances the 2D

Particularly noteworthy is the observation of a critical frequency below which an observed MT tensor becomes sufficiently ill-conditioned for the surficial strike to become evident. The example in Fig. 8 shows this behaviour to apply at frequencies below 10 Hz. That a frequency-independent model produces a frequency-dependent result in this way requires

Remaining problems for which proofs would augment the subject are the evident prohibition of negative *TE* values, the evident non-observation of antiparallel radial arms, and the evident non-observation of MT tensor data with negative determinant values (corresponding to Mohr

The author thanks Peter Milligan, Chris Phillips and John Weaver for many valuable discussions over recent years on the subject matter of this chapter. Data for the NQ142

Bahr, K. (1988). Interpretation of the magnetotelluric impedance tensor: regional induction

Berdichevsky, M. N. & Dmitriev, V. I. (2008). *Models and Methods of Magnetotellurics*, Springer. Bibby, H. M., Caldwell, T. G. & Brown, C. (2005). Determinable and non-determinable

Caldwell, T. G., Bibby, H. M. & Brown, C. (2004). The magnetotelluric phase tensor, *Geophys.*

Eggers, D. E. (1982). An eigenstate formulation of the magnetotelluric impedance tensor,

Groom, R. W. & Bailey, R. C. (1989). Decomposition of magnetotelluric impedance tensors

Gubbins, D. & Herrero-Bervera, E. (eds) (2007). *Encyclopedia of Geomagnetism and*

Hobbs, B. A. (1992). Terminology and symbols for use in studies of electromagnetic induction

LaTorraca, G. A., Madden, T. R. & Korringa, J. (1986). An analysis of the magnetotelluric

parameters of galvanic distortion in magnetotellurics, *Geophys. J. Int.* 163: 915 – 930.

in the presence of local three-dimensional galvanic distortion, *J. Geophys. Res.*

impedance for three-dimensional conductivity structures, *Geophysics* 51: 1819 – 1829.

and local telluric distortion, *J. Geophys.* 62: 119 – 127.

*Paleomagnetism*, Encyclopedia of Earth Sciences, Springer.

in the Earth, *Surv. Geophys.* 13: 489 – 515.

indicate regional geologic strike.

care in interpretation.

**15. Acknowledgements**

**16. References**

regional strike is found in this straightforward way.

circles prohibited from enclosing their origins of axes).

example were provided by Geoscience Australia.

*J. Int.* 158: 457 – 469.

94: 1913 – 1925.

*Geophysics* 47: 1204 – 1214.

By Equations 33, 34, 35 and 36, the values of *θep*, *θhp*, Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* are evaluated as 51.26◦, 24.70◦, 6.36 and −1.88 respectively. As before, by multiplying out the right-hand side, it can be checked that these values do indeed satisfy Equation 32. Also these values may be read off Fig. 12 as 51◦, 25◦, 6.3 and −1.8, but now difficulties arising from the negative determinant are evident. Taking the direction of Υ*<sup>p</sup>* from the origin to define positive, the value of Ψ*<sup>p</sup>* is seen to be negative.

Put into a standard computing routine, the SVD returned for the matrix in Equation 70 is

$$\mathbf{Z}\_p = \begin{bmatrix} .6257 & .7800 \\ .7800 & .6257 \end{bmatrix} \begin{bmatrix} 6.3592 & 0 \\ 0 & 1.8870 \end{bmatrix} \begin{bmatrix} -.4179 & .9085 \\ -.9085 & -.4179 \end{bmatrix} \tag{71}$$

which may also be obtained formally by determining the eigenvalues and eigenvectors of **Z<sup>T</sup>** *<sup>p</sup>***Z***p*. To give it the form of Equation 32, Equation 71 may be expressed

$$\mathbf{Z}\_p = \begin{bmatrix} .6257 & .7800 \\ .7800 & .6257 \end{bmatrix} \begin{bmatrix} 0 & 6.3592 \\ -1.8870 & 0 \end{bmatrix} \begin{bmatrix} .9085 & .4179 \\ -.4179 & .9085 \end{bmatrix} \tag{72}$$

and then written

$$\mathbf{Z}\_p = \begin{bmatrix} \cos 51.26^\circ & \sin 51.26^\circ \\ \sin 51.26^\circ - \cos 51.26^\circ \end{bmatrix} \begin{bmatrix} 0 & 6.3592 \\ -1.8870 & 0 \end{bmatrix} \begin{bmatrix} \cos 24.70^\circ & \sin 24.70^\circ \\ -\sin 24.70^\circ \cos 24.70^\circ \end{bmatrix} \tag{73}$$

The rows of the third matrix [cos 24.70◦, sin 24.70◦; − sin 24.70◦, cos 24.70◦] define unit vectors *v*<sup>1</sup> and *v*<sup>2</sup> for *H*� *<sup>x</sup>*, *H*� *<sup>y</sup>* as shown in Fig. 12. As for Fig. 6, these indicate a rotation of the *Hx*, *Hy* observing axes.

However, the columns of the first matrix [cos 51.26◦, sin 51.26◦; sin 51.26◦, − cos 51.26◦] define unit vectors *u*<sup>1</sup> and *u*<sup>2</sup> also as shown in Fig. 12. These unit vectors do not now represent just a simple rotation of the *Ex*, *Ey* axes, but also the reflection of one of them. The negative value of Ψ*<sup>p</sup>* given by Equations 35 and 36 has been cast by the SVD analysis as a positive singular value in the second matrix in Equation 71. To compensate for this convention, the direction of one of the axes (originating as the direction of an eigenvector of **U** in Equation 29) has been reversed. Any model of distortion rotating the electric fields is thus contravened. An extra physical phenomenon is required to explain the reflection of one of the rotated electric field axes. Pending such an explanation, a negative principal value (−1.8870 in the present example) must be regarded with great caution.

In normal SVD formalism the change of sign of an eigenvector, as demonstrated in the above example, may occur with little comment. In the present case however it has had the profound effect of destroying, after rotation, the "right-handedness" of the observing axes of the MT data.

#### **14. Conclusion**

Basic 2 x 2 matrices arise commonly in linear algebra, and Mohr diagrams have wide application in displaying their properties. Complementing the usual algebraic approach, they may help the student understand especially anomalous data.

Such Mohr diagrams have proved useful in checking MT data for basic errors arising in the manipulation of time-series; in checking hypothetical examples for realistic form; and in checking for errors which are demonstrated by negative determinants.

Invariants of an observed tensor under axes rotation are evident from inspection of such Mohr diagrams. The invariants may be used as indicators of the geologic dimensionality of the tensor; they may also be useful in the inversion and interpretation of MT data. Mohr diagrams of the MT phase tensor also show invariants which measure geologic dimensionality, and indicate regional geologic strike.

Attention has been given to monitoring the singularity of an MT tensor with a condition number, as tensors which approach singularity are common. Carrying out SVD on observed data with high condition numbers can indicate a direction of surficial 2D strike, for the simple case of the surficial distortion of a regional 2D MT response. In some circumstances the 2D regional strike is found in this straightforward way.

Particularly noteworthy is the observation of a critical frequency below which an observed MT tensor becomes sufficiently ill-conditioned for the surficial strike to become evident. The example in Fig. 8 shows this behaviour to apply at frequencies below 10 Hz. That a frequency-independent model produces a frequency-dependent result in this way requires care in interpretation.

Remaining problems for which proofs would augment the subject are the evident prohibition of negative *TE* values, the evident non-observation of antiparallel radial arms, and the evident non-observation of MT tensor data with negative determinant values (corresponding to Mohr circles prohibited from enclosing their origins of axes).

#### **15. Acknowledgements**

The author thanks Peter Milligan, Chris Phillips and John Weaver for many valuable discussions over recent years on the subject matter of this chapter. Data for the NQ142 example were provided by Geoscience Australia.

#### **16. References**

24 Will-be-set-by-IN-TECH

By Equations 33, 34, 35 and 36, the values of *θep*, *θhp*, Υ*<sup>p</sup>* and Ψ*<sup>p</sup>* are evaluated as 51.26◦, 24.70◦, 6.36 and −1.88 respectively. As before, by multiplying out the right-hand side, it can be checked that these values do indeed satisfy Equation 32. Also these values may be read off Fig. 12 as 51◦, 25◦, 6.3 and −1.8, but now difficulties arising from the negative determinant are evident. Taking the direction of Υ*<sup>p</sup>* from the origin to define positive, the value of Ψ*<sup>p</sup>* is

Put into a standard computing routine, the SVD returned for the matrix in Equation 70 is

6.3592 0

which may also be obtained formally by determining the eigenvalues and eigenvectors of

 0 6.3592 −1.8870 0

 0 6.3592 −1.8870 0

The rows of the third matrix [cos 24.70◦, sin 24.70◦; − sin 24.70◦, cos 24.70◦] define unit vectors

However, the columns of the first matrix [cos 51.26◦, sin 51.26◦; sin 51.26◦, − cos 51.26◦] define unit vectors *u*<sup>1</sup> and *u*<sup>2</sup> also as shown in Fig. 12. These unit vectors do not now represent just a simple rotation of the *Ex*, *Ey* axes, but also the reflection of one of them. The negative value of Ψ*<sup>p</sup>* given by Equations 35 and 36 has been cast by the SVD analysis as a positive singular value in the second matrix in Equation 71. To compensate for this convention, the direction of one of the axes (originating as the direction of an eigenvector of **U** in Equation 29) has been reversed. Any model of distortion rotating the electric fields is thus contravened. An extra physical phenomenon is required to explain the reflection of one of the rotated electric field axes. Pending such an explanation, a negative principal value (−1.8870 in the present

In normal SVD formalism the change of sign of an eigenvector, as demonstrated in the above example, may occur with little comment. In the present case however it has had the profound effect of destroying, after rotation, the "right-handedness" of the observing axes of the MT

Basic 2 x 2 matrices arise commonly in linear algebra, and Mohr diagrams have wide application in displaying their properties. Complementing the usual algebraic approach, they

Such Mohr diagrams have proved useful in checking MT data for basic errors arising in the manipulation of time-series; in checking hypothetical examples for realistic form; and

*<sup>y</sup>* as shown in Fig. 12. As for Fig. 6, these indicate a rotation of the *Hx*, *Hy*

0 1.8870

 <sup>−</sup>.4179 .9085 −.9085 −.4179

> .9085 .4179 −.4179 .9085

 cos 24.70◦ sin 24.70◦ − sin 24.70◦ cos 24.70◦ (71)

(72)

(73)

seen to be negative.

and then written

*v*<sup>1</sup> and *v*<sup>2</sup> for *H*�

observing axes.

data.

**14. Conclusion**

**Z***<sup>p</sup>* = 

**Z<sup>T</sup>**

**Z***<sup>p</sup>* = 

**Z***<sup>p</sup>* = 

*<sup>x</sup>*, *H*�

.6257 .7800 .7800 −.6257

.6257 .7800 .7800 −.6257

cos 51.26◦ sin 51.26◦ sin 51.26◦ − cos 51.26◦

example) must be regarded with great caution.

may help the student understand especially anomalous data.

in checking for errors which are demonstrated by negative determinants.

*<sup>p</sup>***Z***p*. To give it the form of Equation 32, Equation 71 may be expressed

Bahr, K. (1988). Interpretation of the magnetotelluric impedance tensor: regional induction and local telluric distortion, *J. Geophys.* 62: 119 – 127.

Berdichevsky, M. N. & Dmitriev, V. I. (2008). *Models and Methods of Magnetotellurics*, Springer.


**5** 

*Italy* 

**Heat and SO2 Emission Rates** 

Letizia Spampinato and Giuseppe Salerno *Istituto Nazionale di Geofisica e Vulcanologia, Osservatorio Etneo, sezione di Catania, Catania,* 

 **of Masaya, Nicaragua** 

 **at Active Volcanoes – The Case Study** 

The necessity of understanding volcanic phenomena so as to assist hazard assessment and risk management, has led to development of a number of techniques for the tracking of volcanic events so as to support forecasting efforts. Since 1980s scientific community has progressively drifted research and surveillance at active volcanoes by integrated approach. Nowadays, volcano observatories over the world record and integrate real or near-real time data for monitoring and understanding volcano behaviour. Among the geophysical, geochemical, and volcanological parameters, the tracking of temperature changes at several volcanic features (e.g. open-vent systems, eruptive vents, fumaroles) and variations in sulphur dioxide flux and concentration at volcanic plumes are key factors for studying and

Temperature is one of the first parameters that have been considered in understanding the nature of volcanoes and their eruptions. Thermal anomalies have proved to be precursors of a number of eruptive events (e.g. Andronico et al., 2005; Dean et al., 2004; Dehn et al., 2002), and once an eruption begins, temperature plays a major role in lava flow emplacement and lava field development (e.g. Ball et al., 2008; Calvari et al., 2010; Lodato et al., 2007). At active volcanoes, temperature has been measured by direct and indirect methodologies (**Fig. 1a, c**). Direct measurements represent the traditional thermal monitoring carried out at fumaroles, hot springs, molten lava bodies, and crater lakes, using thermocouples (e.g. Aiuppa et al., 2006; Corsaro & Miraglia, 2005). Indirect measurements, also known as thermal remote sensing, can be performed by satellite, ground, and airborne surveys (e.g. Calvari et al., 2006; Spampinato et al., 2011; Wright et al., 2010). Owing to the danger of most kinds of eruption, and the need of monitoring inaccessible areas on volcanoes (e.g. Wright & Pilger, 2008), indirect measurements are especially attractive. Among them, thermal imagery is one of the most widespread and results from the capability to detect the infrared radiation emitted from the surface of hot bodies, and to provide the radiometric map of heat distribution of the body's surface (Spampinato et al., 2011). This has been of primary importance for capturing the evolution of thermal anomalies, which shed light on magma movements at shallow depths (e.g. Calvari et al., 2005). While magma is rising, hot gases

**1. Introduction** 

monitoring active volcanoes.

