Chapter 2: The Geomagnetic Field

One of the major efforts in paleomagnetism has been to study ancient geomagnetic fields. Because human measurements extend back about a millennium, measurement of “accidental” records provided by archaeological or geological materials remains the only way to investigate ancient field behavior. Therefore it is useful for students of paleomagnetism to understand something about the present geomagnetic field. In this chapter we review the general properties of the Earth’s magnetic field.

The part of the geomagnetic field of interest to paleomagnetists is generated by convection currents in the liquid outer core of the Earth which is composed of iron, nickel and some unknown lighter component(s). The source of energy for this convection is not known for certain, but is thought to be partly from cooling of the core and partly from the buoyancy of the iron/nickel liquid outer core caused by freezing out of the pure iron inner core. Motions of this conducting fluid are controlled by the buoyancy of the liquid, the spin of the Earth about its axis and by the interaction of the conducting fluid with the magnetic field (in a horribly non-linear fashion). Solving the equations for the fluid motions and resulting magnetic fields is a challenging computational task. Recent numerical models, however, show that such magnetohydrodynamical systems can produce self-sustaining dynamos which create enormous external magnetic fields.

2.1Components of magnetic vectors¶

The magnetic field of a dipole aligned along the spin axis and centered in the Earth (a so-called geocentric axial dipole, or GAD) is shown in Figure 2.1a. [See Chapter 1 for a derivation of how to find the radial and tangential components of such a field.] By convention, the sign of the Earth’s dipole is negative, pointing toward the south pole as shown in Figure 2.1a and magnetic field lines point toward the north pole. They point downward in the northern hemisphere and upward in the southern hemisphere.

Although dominantly dipolar, the geomagnetic field is not perfectly modeled by a geocentric axial dipole, but is somewhat more complicated (see Figure 2.1b). At the point on the surface labeled ‘P’, the geomagnetic field points nearly north and down at an angle of approximately 60°. Vectors in three dimensions are described by three numbers and in many paleomagnetic applications, these are two angles ($D$ and $I$) and the strength ($B$) as shown in Figure 2.1b and c. The angle from the horizontal plane is the inclination $I$; it is positive downward and ranges from +90° for straight down to -90° for straight up. If the geomagnetic field were that of a perfect GAD field, the horizontal component of the magnetic field ($B_H$ in Figure 2.1b) would point directly toward geographic north. In most places on the Earth there is a deflection away from geographic north and the angle between geographic and magnetic north is the declination, $D$ (see Figure 2.1c). $D$ is measured positive clockwise from North and ranges from $0 \rightarrow 360^\circ$. [Westward declinations can also be expressed as negative numbers, i.e., 350° = -10°.] The vertical component ($B_V$ in Figure 2.1b,c) of the geomagnetic field at P, is given by

and the horizontal component $B_H$ (Figure 2.1c) by

$B_H$ can be further resolved into north and east components ($B_N$ and $B_E$ in Figure 2.1c) by

Depending on the particular problem, some coordinate systems are more suitable to use because they have the symmetry of the problem built into them. We have just defined a coordinate system using two angles and a length ($B,D,I$) and the equivalent Cartesian coordinates of ($B_N, B_E, B_V$). We will need to convert among them at will. There are many names for the Cartesian coordinates. In addition to north, east and down, they could also be $x,y,z$ or even $x_1, x_2$ and $x_3$. The convention used in this book is that axes are denoted $\mathbf{X}_1, \mathbf{X}_2, \mathbf{X}_3$, while the components along the axes are frequently designated $x_1, x_2, x_3$. In the geographic frame of reference, positive $\mathbf{X}_1$ is to the north, $\mathbf{X}_2$ is east and $\mathbf{X}_3$ is vertically down in keeping with the right-hand rule. To convert from Cartesian coordinates to angular coordinates ($B,D,I$):

Be careful of the sign ambiguity of the tangent function. You may well end up in the wrong quadrant and have to add 180°; this will happen if both $x_1$ and $x_2$ are negative. In most computer languages, there is a function atan2 which takes care of this, but most hand calculators will not. Remember that most computer languages expect angles to be given in radians, not degrees, so multiply degrees by $\pi/180$ to convert to radians. Note also that in place of $\mathbf{B}$ for magnetic induction with units of tesla as a measure of vector length, (see Chapter 1), we could also use $\mathbf{H}$, $\mathbf{M}$ (both Am$^{-1}$) or $\mathbf{m}$ (Am$^2$) for magnetic field, magnetization or magnetic moment respectively.

2.2Reference magnetic field¶

We can measure declination, inclination and intensity at different places around the globe, but not everywhere all the time. Yet it is often handy to be able to predict what these components are. For example, it is extremely useful to know what the deviation is between true North and declination in order to find our way with maps and compasses. In principle, magnetic field vectors can be derived from the magnetic potential $\psi_m$ as we showed in Chapter 1. For an axial dipolar field, there is but one scalar coefficient (the magnetic moment $\mathbf{m}$ of a dipole source). For the geomagnetic field, there are many more coefficients, including not just an axial dipole aligned with the spin axis, but two orthogonal equatorial dipoles and a whole host of more complicated sources such as quadrupoles, octupoles and so on. A list of coefficients associated with these sources allows us to calculate the magnetic field vector anywhere outside of the source region. In this section, we outline how this might be done.

As we learned in Chapter 1, the magnetic field at the Earth’s surface can be calculated from the gradient of a scalar potential field ($\mathbf{H}=-\nabla \psi_m$), and this scalar potential field satisfies Laplace’s Equation:

For the geomagnetic field (ignoring external sources of the magnetic field which are in any case small and transient), the potential equation can be written as:

where $a$ is the radius of the Earth ($6.371 \times 10^6$ m). In addition to the radial distance $r$ and the angle away from the pole $\theta$, there is $\phi$, the angle around the equator from some reference, say, the Greenwich meridian. Here, $\theta$ is the co-latitude and $\phi$ is the longitude. The $g_l^m$s and $h_l^m$s are the Gauss coefficients (degree $l$ and order $m$) for hypothetical sources at radii less than $a$ calculated for a particular year. These are normally given in units of nT. The $P_l^m$s are wiggly functions called partially normalized Schmidt polynomials of the argument $\cos \theta$. These are closely related to the associated Legendre polynomials. [When $m=0$ the Schmidt and Legendre polynomials are identical.] The first few $P_l^m$s are:

and are shown in Figure 2.2.

To get an idea of how the Gauss coefficients in the potential relate to the associated magnetic fields, we show three examples in Figure 2.3. We plot the inclinations of the vector fields that would be produced by the terms with $g_1^0, g_2^0$ and $g_3^0$ respectively. These are the axial ($m=0$) dipole ($l=1$), quadrupole ($l=2$) and octupole ($l=3$) terms. The associated potentials for each harmonic are shown in the insets.

In general, terms for which the difference between the subscript ($l$) and the superscript ($m$) is odd (e.g., the axial dipole $g_1^0$ and octupole $g_3^0$) produce magnetic fields that are antisymmetric about the equator, while those for which the difference is even (e.g., the axial quadrupole $g_2^0$) have symmetric fields. In Figure 2.3a we show the inclinations produced by a purely dipolar field of the same sign as the present day field. The inclinations are all positive (down) in the northern hemisphere and negative (up) in the southern hemisphere. In contrast, inclinations produced by a purely quadrupolar field (Figure 2.3b) are down at the poles and up at the equator. The map of inclinations produced by a purely axial octupolar field (Figure 2.3c) are again asymmetric about the equator with vertical directions of opposite signs at the poles separated by bands with the opposite sign at mid-latitudes.

As noted before, there is not one, but three dipole terms in Equation 2.6, the axial term ($g_1^0$) and two equatorial terms ($g_1^1$ and $h_1^1$). Therefore, the total dipole contribution is the vector sum of these three or $\sqrt{{g_1^0}^2+{g_1^1}^2+{h_1^1}^2}$. The total quadrupole contribution ($l=2$) combines five coefficients and the total octupole ($l=3$) contribution combines seven coefficients.

So how do we get this marvelous list of Gauss coefficients? If you want to know the details, please refer to Langel (1987). We will just give a brief introduction here. Recalling Chapter 1, once the scalar potential $\psi_m$ is known, the components of the magnetic field can be calculated from it. We solved this for the radial and tangential field components ($H_r$ and $H_{\theta}$) in Chapter 1. We will now change coordinate and unit systems and introduce a third dimension (because the field is not perfectly dipolar). The north, east, and vertically down components are related to the potential $\psi_m$ by:

where $r$, $\theta$, $\phi$ are radius, co-latitude (degrees away from the North pole) and longitude, respectively. Here, $B_V$ is positive down, $B_E$ is positive east, and $B_N$ is positive to the north, the opposite of $H_r$ and $H_{\theta}$ as defined in Chapter 1. Note that Equation 2.8 is in units of induction, not Am$^{-1}$ if the units for the Gauss coefficients are in nT, as is the current practice.

Going backwards, the Gauss coefficients are determined by fitting Equations 2.8 and 2.6 to observations of the magnetic field made by magnetic observatories or satellite for a particular time. The International (or Definitive) Geomagnetic Reference Field or I(D)GRF, for a given time interval is an agreed upon set of values for a number of Gauss coefficients and their time derivatives. IGRF (or DGRF) models and programs for calculating various components of the magnetic field are available on the internet from the National Geophysical Data Center; the address is http://www.ngdc.noaa.gov. There is also a program igrf.py included in the PmagPy package (see igrf.py documentation).

In practice, the Gauss coefficients for a particular reference field are estimated by least-squares fitting of observations of the geomagnetic field. You need a minimum of 48 observations to estimate the coefficients to $l=6$. Nowadays, we have satellites which give us thousands of measurements and the list of generation 10 of the IGRF for 2005 goes to $l=13$.

In order to get a feel for the importance of the various Gauss coefficients, take a look at Table 2.1, which has the Schmidt quasi-normalized Gauss coefficients for the first six degrees from the IGRF for 2005. The power at each degree is the average squared field per spherical harmonic degree over the Earth’s surface and is calculated by $R_l=\sum_m (l+1)[(g_l^m)^2 + (h_l^m)^2]$ Lowes, 1974. The so-called Lowes spectrum is shown in Figure 2.4. It is clear that the lowest order terms (degree one) totally dominate, constituting some 90% of the field. This is why the geomagnetic field is often assumed to be equivalent to a magnetic field created by a simple dipole at the center of the Earth.

2.3Geocentric axial dipole (GAD) and other poles¶

The beauty of using the geomagnetic potential field is that the vector field can be evaluated anywhere outside the source region. Using the values for a given reference field in Equations 2.6 and 2.8, we can calculate values of $B, D$ and $I$ at any location on Earth. Figure 2.1b shows the lines of flux predicted from the 2005 IGRF from the core-mantle boundary up. We can see that the field becomes simpler and more dipolar as we move from the core mantle boundary to the surface. Yet, there is still significant non-dipolar structure in the geomagnetic field even at the Earth’s surface.

We can recast the vectors at the surface of the Earth into maps of components as shown in Figure 2.5a,b. We show the potential in Figure 2.5c for comparison with that of a pure dipole (inset to Figure 2.3a). These maps illustrate the fact that the field is a complicated function of position on the surface of the Earth. The intensity values in Figure 2.5a are, in general, highest near the poles (~60 $\mu$T) and lowest near the equator (~30 $\mu$T), but the contours are not straight lines parallel to latitude as they would be for a field generated strictly by a geocentric axial dipole (GAD) (e.g, Figure 2.1a). Similarly, a GAD would produce lines of inclination that vary in a regular way from -90° to +90° at the poles, with 0° at the equator; the contours would parallel the lines of latitude. Although the general trend in inclination shown in Figure 2.5b is similar to the GAD model, the field lines are more complicated, which again suggests that the field is not perfectly described by a geocentric bar magnet.

Perhaps the most important result of spherical harmonic analysis for our purposes is that the field at the Earth’s surface is dominated by the degree one terms ($l=1$) and the external contributions are very small. The first order terms can be thought of as geocentric dipoles that are aligned with three different axes: the spin axis ($g_1^0$) and two equatorial axes that intersect the equator at the Greenwich meridian ($h_1^0$) and at 90° East ($h_1^1$). The vector sum of these geocentric dipoles is a dipole that is currently inclined by about 10° to the spin axis. The axis of this best-fitting dipole pierces the surface of the Earth at the circle in Figure 2.6. This point and its antipode are called geomagnetic poles. Points at which the field is vertical ($I = \pm 90^{\circ}$ shown by a square in Figure 2.6) are called magnetic poles, or sometimes dip poles. These poles are distinguishable from the geographic poles, where the spin axis of the Earth intersects its surface. The northern geographic pole is shown by a star in Figure 2.6.

It turns out that when averaged over sufficient time, the geomagnetic field actually does seem to be approximately a GAD field, perhaps with a pinch of $g_2^0$ thrown in (see e.g., Merrill et al. (1996)). The GAD model of the field will serve as a useful crutch throughout our discussions of paleomagnetic data and applications. Averaging ancient magnetic poles over enough time to average out secular variation (thought to be 10$^4$ or 10$^5$ years) gives what is known as a paleomagnetic pole; this is usually assumed to be co-axial with the Earth’s geographic pole (the spin axis).

Because the geomagnetic field is axially dipolar to a first approximation, we can write:

Note that $g_1^0$ is given in nT in Table 2.1. Thus, from Equation 2.9,

Given some latitude $\lambda$ on the surface of the Earth in Figure 2.1a and using the equations for $B_V$ and $B_N$, we find that:

This equation is sometimes called the dipole formula and shows that the inclination of the magnetic field is directly related to the co-latitude ($\theta$) for a field produced by a geocentric axial dipole (or $g_1^0$). The dipole formula allows us to calculate the latitude of the measuring position from the inclination of the (GAD) magnetic field, a result that is fundamental in plate tectonic reconstructions. The intensity of a dipolar magnetic field is also related to (co)latitude because:

The dipole field intensity has changed by more than an order of magnitude in the past and the dipole relationship of intensity to latitude turns out to be not useful for tectonic reconstructions.

2.4Plotting magnetic directional data¶

Magnetic field and magnetization directions can be visualized as unit vectors anchored at the center of a unit sphere. Such a unit sphere is difficult to represent on a 2-D page. There are several popular projections, including the Lambert equal area projection which we will be making extensive use of in later chapters. The principles of construction of the equal area projection are covered in Section B.1.

In general, regions of equal area on the sphere project as equal area regions on this projection, as the name implies. Plotting directional data in this way enables rapid assessment of data scatter. A drawback of this projection is that circles on the surface of a sphere project as ellipses. Also, because we have projected a vector onto a unit sphere, we have lost information concerning the magnitude of the vector. Finally, lower and upper hemisphere projections must be distinguished with different symbols. The paleomagnetic convention is: lower hemisphere projections (downward directions) use solid symbols, while upper hemisphere projections are open.

The dipole formula allows us to convert a given measurement of $I$ to an equivalent magnetic co-latitude $\theta_m$:

If the field were a simple GAD field, $\theta_m$ would be a reasonable estimate of $\theta$, but non-GAD terms can invalidate this assumption. To get a feel for the effect of these non-GAD terms, we consider first what would happen if we took random measurements of the Earth’s present field (see Figure 2.7). We evaluated the directions of the magnetic field using the IGRF for 2005 at 200 positions on the globe (shown in Figure 2.7a). These directions are plotted in Figure 2.7b using the paleomagnetic convention of open symbols pointing up and closed symbols pointing down. In Figure 2.7c, we plot the inclinations as a function of latitude. As expected from a predominantly dipolar field, inclinations cluster around the values for a geocentric axial dipolar field but there is considerable scatter and interestingly the scatter is larger in the southern hemisphere than in the northern one. This is related to the low intensities beneath South America and the Atlantic region seen in Figure 2.5a.

2.4.1$D', I'$ transformation¶

Often we wish to compare directions from distant parts of the globe. There is an inherent difficulty in doing so because of the large variability in inclination with latitude. In such cases it is appropriate to consider the data relative to the expected direction (from GAD) at each sampling site. For this purpose, it is useful to use a transformation whereby each direction is rotated such that the direction expected from a geocentric axial dipole field (GAD) at the sampling site is the center of the equal area projection. This is accomplished as follows:

Each direction is converted to Cartesian coordinates ($x_{i}$) by:

These are rotated to the new coordinate system ($x'_{i}$, see Section A.3.5.2) by:

where $I_{d}$ = the inclination expected from a GAD field ($\tan I_{d}=2\tan \lambda$), $\lambda$ is the site latitude, and $\alpha$ is the inclination of the paleofield vector projected onto the N-S plane ($\alpha = \tan^{-1} (x_{3}/x_{1})$). The $x'_{i}$ are then converted to $D',I'$ by Equation 2.4.

In Figure 2.8a we show the geomagnetic field vectors evaluated at random longitudes along a latitude band of 45°N. The vectors are shown in their Cartesian coordinates of North, East and Down. In Figure 2.8b we show what happens when we rotate the coordinate system to peer down the direction expected from an axial dipolar field at 45°N (which has an inclination of 63°). The vectors circle about the expected direction. Finally, we see what happens to the directions shown in Figure 2.7b after the $D',I'$ transformation in Figure 2.8. These are unit vectors projected along the expected direction for each observation in Figure 2.7a. Comparing the equal area projection of the directions themselves (Figure 2.7b) to the transformed directions (Figure 2.8c), we see that the latitudinal dependence of the inclinations has been removed.

2.4.2Virtual geomagnetic poles¶

We are often interested in whether the geomagnetic pole has changed, or whether a particular piece of crust has rotated with respect to the geomagnetic pole. Yet, what we observe at a particular location is the local direction of the field vector. Thus, we need a way to transform an observed direction into the equivalent geomagnetic pole.

In order to remove the dependence of direction merely on position on the globe, we imagine a geocentric dipole which would give rise to the observed magnetic field direction at a given latitude ($\lambda$) and longitude ($\phi$). The virtual geomagnetic pole (VGP) is the point on the globe that corresponds to the geomagnetic pole of this imaginary dipole (Figure 2.9a).

Paleomagnetists use the following conventions: $\phi$ is measured positive eastward from the Greenwich meridian and ranges from $0 \rightarrow 360^\circ$; $\theta$ is measured from the North pole and goes from $0 \rightarrow 180^\circ$. Of course $\theta$ relates to latitude, $\lambda$ by $\theta = 90 - \lambda$. $\theta_m$ is the magnetic co-latitude and is given by Equation 2.13. Be sure not to confuse latitudes and co-latitudes. Also, be careful with declination. Declinations between 180° and 360° are equivalent to $D$ - 360 ° which are counter-clockwise with respect to North.

The first step in the problem of calculating a VGP is to determine the magnetic co-latitude $\theta_m$ by Equation 2.13 which is defined in the dipole formula (Equation 2.13). The declination $D$ is the angle from the geographic North Pole to the great circle joining the observation site $S$ and the pole $P$, and $\Delta \phi$ is the difference in longitudes between P and S, $\phi_p-\phi_s$. Now we use some tricks from spherical trigonometry as reviewed in Section A.3.1.

We can locate VGPs using the law of sines and the law of cosines. The declination $D$ is the angle from the geographic North Pole to the great circle joining $S$ and $P$ (see Figure 2.9) so:

which allows us to calculate the VGP co-latitude $\theta_p$. The VGP latitude is given by:

so $90>\lambda_p>0$ in the northern hemisphere and $0<\lambda_p<90$ in the southern hemisphere.

To determine $\phi_p$, we first calculate the angular difference between the pole and site longitude $\Delta \phi$.

If $\cos \theta_m \geq \cos \theta_s \cos \theta_p$, then $\phi_p = \phi_s+\Delta \phi$. However, if $\cos \theta_m < \cos \theta_s \cos \theta_p$ then $\phi_p=\phi_s + 180 -\Delta \phi$.

Now we can convert the directions in Figure 2.7b to VGPs (see Figure 2.9c). The grouping of points is much tighter in Figure 2.9c than in the equal area projection because the effect of latitude variations in dipole fields has been removed. If a number of VGPs are averaged together, the average pole position is called a “paleomagnetic pole”. How to average poles and directions is the subject of Chapter 11 and Chapter 12.

The procedure for calculating a direction from a VGP is a similar procedure to that for calculating the VGP from the direction. Magnetic colatitude $\theta_m$ is calculated in exactly the same way as before and yields inclination from the dipole formula. The declination can be calculated by solving for $D$ in Equation 2.16 as:

This equation works most of the time, but breaks down under some circumstances, for example, when the pole latitude is further to the south than the site latitude. The following algorithm works in the more general case:

where $C = |1- (\cos D)^2|$. Also, if $-90 < \Delta \phi <0$ or if $\Delta \phi > 180$, then $D = 360 - D$.

2.4.3Virtual dipole moment¶

As pointed out earlier, magnetic intensity varies over the globe in a similar manner to inclination. It is often convenient to express paleointensity values in terms of the equivalent geocentric dipole moment that would have produced the observed intensity at a specific (paleo)latitude. Such an equivalent moment is called the virtual dipole moment (VDM) by analogy to the VGP (see Figure 2.9a). First, the magnetic (paleo)co-latitude $\theta_m$ is calculated as before from the observed inclination and the dipole formula of Equation 2.11. Then, following the derivation of Equation 2.12, we have

Sometimes the site co-latitude as opposed to magnetic co-latitude is used in the above equation, giving a virtual axial dipole moment (VADM; see Figure 2.9d).

SUPPLEMENTAL READINGS: Merrill et al. (1996), Chapters 1 & 2

2.5Problems¶

For this and future problem sets, you will need the PmagPy package (see section in the Preface at the beginning of the book). After you have installed this and properly set your path, you can import the functions from PmagPy using these commands:

import pmagpy.pmag as pmag
import pmagpy.ipmag as ipmag
import pmagpy.pmagplotlib as pmagplotlib

Please consult the Jupyter notebook PmagPy.ipynb for more help on using PmagPy functions within a notebook.

Problem 1

a) Write a python script in a Jupyter notebook that converts declination, inclination and intensity to North, East, and Down. Read in the data in the file Chapter_2/ps2_prob1_data.txt. For this the loadtxt function in the Numpy module will come in handy.

b) Choose 10 random spots on the surface of the earth. You can use the pmag.get_unf to generate a list for you. Then use the ipmag.igrf function to evaluate the declination, inclination and intensity at each of these locations in January 2006. As with all PmagPy programs, and functions, you can find out what they do by printing out the doc string: you can find out what they do by getting the help message:

help(ipmag.igrf)

Calls like these generates help messages which will help you to call the function properly.

c) Take the vectors from the output of Problem 1b and convert them to cartesian coordinates, using the script you wrote in Problem 1a.

Problem 2

a) Plot the IGRF directions from Problem 1b on an equal area projection by hand. Use the equal area net provided in Section B.1. Remember that the outer rim is horizontal and the center of the diagram is vertical. Azimuth goes around the rim with clockwise being positive. Put a thumbtack through the equal area (Schmidt) net and place a piece of tracing paper on the thumbtack. Mark the top of the stereonet with a tick mark on the tracing paper.

To plot a direction, rotate the tick mark of the tracing paper around counter clockwise until the top of the paper is rotated by the declination of the direction. Then count tick marks toward the center from the outer rim (the horizontal) to the inclination angle, plot the point, and rotate back so that the tick is North again. Put all your points on the diagram.

b) Now use the ipmag functions plot_net and plot_di. or write your own! Both plots should look the same....

Problem 3

You went to Wyoming (112° W and 36° N) to sample some Cretaceous rocks. You measured a direction with a declination of 345° and an inclination of 47°.

a) What direction would you expect from the present (GAD) field?

b) What is the virtual geomagnetic pole position corresponding to the direction you actually measured? [Hint: Use the function pmag.dia_vgp in the PmagPy module or for a challenge, write your own!]