Appendix A Historical geomagnetic field model gufm1

49 Pages·2012·3.78 MB·English

249
Appendix A
Historical geomagnetic (cid:12)eld model gufm1
A.1 Historical observations of Earth’s magnetic (cid:12)eld
Direct observations of the geomagnetic (cid:12)eld with the extensive geographical coverage
required for global (cid:12)eld modelling exist from the 16th century onwards from a variety of
sources:
(i) Maritime observations
(ii) Observatory measurements (after 1830)
(iii) Magnetic surveys (18th, 19th and 20th century)
(iv) Satellite magnetometer data (POGO late 1960s, Magsat 1980)
Jackson et al. (2000) constructed a time-dependent (cid:12)eld model known as gufm1 based
on a new compilation of all suitable historical observations. Model gufm1 represents a
signi(cid:12)cant advance on previous historical models such as the ufm model of Bloxham and
Jackson (1992), particularly in the number and spatial distribution of observations in
the 18th and 19th centuries, though there are still changes in coverage with time (see
(cid:12)gure A.1). The temporal change in the number of observations in each (cid:12)ve year period
between 1590 and 1990 is displayed in (cid:12)gure A.2.
A.2 Field modelling methodology
An outline of the inversion techniques used to construct the (cid:12)eld model gufm1 from
historical observations is given here to aid understandingof the limitations of the model.
TheconclusionsfromtheanalysisinChapter3ofthepropertiesof(cid:12)eldfeaturesingufm1
are only trusted if the (cid:12)eld features seen in the model can be considered robust.
Assuming the mantle to be an insulator, then the magnetic (cid:12)eld outside the core can be
expressed as the gradient of a scalar potential V that satis(cid:12)es Laplace’s equation and is
250 Appendix A | gufm1
(a) Locations of declination observations from 1650-1699
(b) Locations of declination observations from 1750-1799
(c) Locations of all observations from 1850-1899
Figure A.1: Changes in spatial distribution of observations in gufm1.
Stars show the location of observations for the (cid:12)fty year periods 1650-1699 in (a), 1750-
1799 in (b), and 1850-1899 in (c). For further details see Jackson et al. (2000) and
Jonkers et al. (2003)
251 Appendix A | gufm1
25000
20000
y
c
n15000
e
u
q
e
r10000
F
5000
0
1600 1700 1800 1900 2000
Date
Figure A.2: Temporal distribution of gufm1 observations.
The number of measurements used in constructing gufm1 in (cid:12)ve year intervals| for
further details see Jackson et al. (2000) and Jonkers et al. (2003)
also a function of time
B = V; (A.1)
(cid:0)r
where 2V = 0; (A.2)
r
1
and V as r ; (A.3)
(cid:0)! r2 (cid:0)!1
where it has also been assumed that there are no (cid:12)eld sources at in(cid:12)nity. Choosing a
spherical polar co-ordinate system (r;(cid:18);(cid:30)) the method of separation of variables is used
to(cid:12)ndasolutiontoLaplace’s equation subjecttotheseboundaryconditions, intheform
of a spherical harmonic expansion
1 l
a l+1
V(t) = a gm(t)cos(m(cid:30))+hm(t)sin(m(cid:30)) Pm(cos(cid:18)); (A.4)
r f l l g l
Xl=1mX=0(cid:16) (cid:17)
whereonlythesolutionscorrespondingtointernal(cid:12)eldsources (foundinmodellingtobe
the dominant terms) have been retained, a=6371.2 km is Earth’s mean radius and Pm
l
are the Schmidt quasi-normalised associated Legendre functions of degree l and order m
with normalisation de(cid:12)ned as
2(cid:25) (cid:25)
Pm(cos(cid:18))Pm0(cos(cid:18)) cosm(cid:30)cosm0(cid:30) sin(cid:18)d(cid:18)d(cid:30)= 24l+(cid:25)1 ifl=l0;m=m0 : (A.5)
l l0 sinm(cid:30)sinm0(cid:30) 0 otherwise
Z Z (cid:26) (cid:27) (cid:26) (cid:27)
0 0
The constants gm(t), hm(t) are expanded in fourth order B-spline basis functions M (t)
l l n
erected on a series of knots such that
gm(t) = gmnM (t); (A.6)
l l n
n
X
252 Appendix A | gufm1
hm(t) = hmnM (t); (A.7)
l l n
n
X
whereM (t)> 0 if t (cid:15) t , t andis zero otherwise.Taking thegradient of equation
n n n+4
j j
((A.4)) we obtain the following expressions for the spherical polar components of the
magnetic (cid:12)eld B: the radial magnetic (cid:12)eld B , the southward magnetic (cid:12)eld B and the
r (cid:18)
eastward magnetic (cid:12)eld B .
(cid:30)
1 l
a l+2
B = (l+1) gm(t)cos(m(cid:30))+hm(t)sin(m(cid:30)) Pm(cos(cid:18)) (A.8)
r r l l l
Xl=1mX=0 (cid:16) (cid:17) (cid:8) (cid:9)
1 l a l+2 dPm(cos(cid:18))
B = gm(t)cos(m(cid:30))+hm(t)sin(m(cid:30)) l (A.9)
(cid:18) (cid:0) r l l d(cid:18)
Xl=1mX=0(cid:16) (cid:17) (cid:8) (cid:9)
1 l
1 a l+2
B = m gm(t)sin(m(cid:30)) hm(t)cos(m(cid:30)) Pm(cos(cid:18))(A.10)
(cid:30) sin(cid:18) r l (cid:0) l l
Xl=1mX=0 (cid:16) (cid:17) (cid:8) (cid:9)
The inverse problem of geomagnetic (cid:12)eld modelling now involves (cid:12)nding the parameters
gmn;hmn that determine the most reasonable (in a sense de(cid:12)ned below) magnetic (cid:12)eld
l l
at the core surface given geomagnetic observations and a priori information on the (cid:12)eld
structure. The prior information includes the fact that the mantle is approximately an
insulator, that the internal (cid:12)eld must originate in the core, the known core radius, and
an upper limit on Ohmic heating from the magnetic (cid:12)eld derived from the induction
equation (Gubbins, 1975).
Spherical harmonic expansions were truncated at degree l = L where L = 14 to obtain a
(cid:12)nite representation of an otherwise in(cid:12)nite dimensional inverse problem. However, the
favoured(cid:12)eldmodelsconvergeatl < Lthankstotheregularisationappliedthatpenalises
spatially complex (cid:12)eld structures (see further discussion of regularisation below). 163
B-splines were employed to represent the time period 1590 to 1990, erected on N = 167
knots (spline end points) placed every 2.5 years. The total number of free parameters in
the model is thus P = N(L(L+2)) = 36512.
Thedataintheinverseproblemweremagnetic(cid:12)eldobservations ofvarioustypes: north-
ward (cid:12)eld component (X), eastward (cid:12)eld component (Y), downward (cid:12)eld component
(Z), declination (D), inclination (I), horizontal (cid:12)eld component (H) and total (cid:12)eld (F).
These are related to the spherical polar (cid:12)eld components (B ;B ;B ) and hence the
r (cid:18) (cid:30)
model parameters by the relations
X = B ; (A.11)
(cid:18)
(cid:0)
Y = B ; (A.12)
(cid:30)
Z = B ; (A.13)
r
(cid:0)
1
H = B2+B2 2 ; (A.14)
(cid:18) (cid:30)
(cid:0) (cid:1)
253 Appendix A | gufm1
1
F = B2+B2 +B2 2 ; (A.15)
(cid:18) (cid:30) r
(cid:0) (cid:1)
Br (cid:25) (cid:25)
I = Arctan (cid:0) where I ; (A.16)
2 13 (cid:0) 2 (cid:20) (cid:20) 2
B2+B2 2
6 (cid:18) (cid:30) 7
4(cid:16) (cid:17) 5
B
(cid:30)
D = Arctan where (cid:25) D (cid:25); (A.17)
B (cid:0) (cid:20) (cid:20)
(cid:20)(cid:0) (cid:18)(cid:21)
Equations((A.11))to((A.13))arelinearinthemodelcoe(cid:14)cients gmn;hmn whilethose
f l l g
in equations ((A.14)) to ((A.17)) are non-linear. This means the whole inverse problem
must be solved iteratively if all the data types are to be included. If the observational
data X;Y;Z;H;F;I;D are listed in a vector d of observations ( = 365 694 in gufm1),
D D
and the model coe(cid:14)cients are listed in a vector m = (g00;g10;h10::) then the forward
1 1 1
problem can be written in matrix form as
d = f(m)+e; (A.18)
where f is the non-linear functional relating the data to the model and e is an error
vector of length which gives the error associated with each observation. Errors vary
D
depending on the method of data acquisition.
The inverse problem of (cid:12)nding the model parameters (spherical harmonic coe(cid:14)cients)
that best represent the core surface (cid:12)eld, given geomagnetic measurements at Earth’s
surface is non-unique. The non-uniqueness arises due to the (cid:12)nite number of data (mag-
netic (cid:12)eld observations at speci(cid:12)c locations and times) being used to (cid:12)nd a continuous
function (the global, time-dependent magnetic (cid:12)eld at the core surface). Geophysically
plausible solutions can be found by using the method of regularised least squares in-
version, that looks for solutions that both (cid:12)t the data as well as possible in a least
squares sense and minimise suitable model norms. Two model norms are employed in
the construction of gufm1, one measuring roughness in the spatial domain and the other
roughness in the temporal domain. The spatial norm is a quadratic norm that is physi-
cally associated with minimising Ohmic dissipation (Gubbins, 1975) and can be written
as
1 te
(cid:9) = (B )dt = mTS(cid:0)1m; (A.19)
r
t t F
e(cid:0) s Zts
L l
(l+1)(2l+1)(2l+3) a 2l+4
where (B ) =4(cid:25) (gm)2+(hm)2 ; (A.20)
F r l c l l
Xl=1 (cid:16) (cid:17) mX=0(cid:2) (cid:3)
The temporal norm is
1 te
(cid:8) = (@2B )2d(cid:10)dt= mTT(cid:0)1m; (A.21)
t t t r
e(cid:0) s Zts ICMB
Measurements(duetocrustal(cid:12)elds,external(cid:12)elds,(cid:12)eldmeasurementerrors,mislocation
errors) were estimated and assumed to be Gaussian (though this may not be totally
254 Appendix A | gufm1
correct | see Walker and Jackson (2001) for a discussion), allowing the construction of
the data covariance matrix C . The model estimate is then that which minimises the
e
objective function (cid:2)
(cid:2)(m)= [d f(m)]TC(cid:0)1[d f(m)]+mTC(cid:0)1m; (A.22)
(cid:0) e (cid:0) m
with the model covariance matrix being
C(cid:0)1 = (cid:21) S(cid:0)1+(cid:21) T(cid:0)1 (A.23)
m S T
(cid:0) (cid:1)
where (cid:21) is the spatial damping parameter (chosen to be 1 10(cid:0)12nT(cid:0)2) and (cid:21) is the
S T
(cid:2)
temporal damping parameter (chosen to be 5 10(cid:0)4nT(cid:0)2yr(cid:0)4). The optimal solution is
(cid:2)
sought by iterating using the scheme in (A.24) below, until convergence on a model with
minimum (cid:2). The scheme is derived by di(cid:11)erentiating (A.22) with respect to m , the
i
ith model parameter, and using a (cid:12)rst order approximation f(m+(cid:14)m) =f(m)+A(cid:14)m
where A = @fi and (cid:14)m is the di(cid:11)erence between the current approximation and the
ij @mj
minimising model
m = m +(ATC(cid:0)1A+C(cid:0)1)(cid:0)1[ATC(cid:0)1(d f(m )) C(cid:0)1m ]: (A.24)
i+1 i e m e (cid:0) i (cid:0) m i
The regularisation ensures that expansions converge before the truncation levels are
reached, so the solutions are insensitive to truncation. The procedure e(cid:11)ectively damps
out most power at degrees 12 and higher. The choice of damping parameters is how-
ever rather subjective|large di(cid:11)erences in model complexity can be produced by vary-
ing them. The values used here were chosen to produce qualitatively similar maps
to those of other authors at (cid:12)xed epochs (refer e.g. Shure et al. (1985)) yielding
a conservative picture of the (cid:12)eld structures required to produce a normalised mis(cid:12)t
(cid:31) = [d(cid:0)f(m)]TC(cid:0)e1[d(cid:0)f(m)] , where D is the total number of data observations, as
D
close qto 1.0 as possible (i.e. (cid:12)tting observations to within estimated errors). The pre-
ferred model gufm1 actually has a normalised mis(cid:12)t of 1.16. More complex models with
more power at higher spherical harmonic degree and larger time variations that (cid:12)t ob-
servations better can be produced, but analysing (cid:12)eld features from these models and
basing interpretations on them increases the chances that the data are being over-(cid:12)t or
e(cid:11)ects due to non-core (cid:12)elds erroneously included.
A.3 Comparison to observatory data
Perhapsthemostrigoroustestofthereliabilityofthemodelgufm1ishowwellitmanages
to (cid:12)tmain (cid:12)eld andsecular variation data recorded at observatories. Here the main (cid:12)eld
and secular variation from 5 observatories are plotted in (cid:12)gure (A.3) and (cid:12)gure (A.4)
respectively. The main (cid:12)eld component typically has little scatter and gufm1 excellently
(cid:12)ts the signal at all the observatories considered. There is considerably more scatter in
255 Appendix A | gufm1
(a) Annual means at Fuquene (MF) (b) Annual means at Mauritius (MF)
(c) Annual means at Niemegk (MF) (d) Annual means at Tuscon (MF)
Figure A.3: Comparison of gufm1 and observatory annualmeans (main(cid:12)eld).
Comparison of the main (cid:12)eld predicted by the model gufm1 and annual means of the
main (cid:12)eld components (X = B , Y = B and Z = B ) measured at observatories in
(cid:18) (cid:30) r
(cid:0) (cid:0)
di(cid:11)erent locations (Fuquene, Mauritius, Niemegk and Tuscon). Observatory data were
obtained from the world data centre. The black squares, triangles and crosses are the
observatory annual means and the red line is the prediction of gufm1. There was a site
change in Mauritius in 1920, observatory results from both before and after the site
change are presented in (b).
256 Appendix A | gufm1
(a) Annual means at Fuquene (SV) (b) Annual means at Mauritius (SV)
(c) Annual means at Niemegk (SV) (d) Annual means at Tuscon (SV)
Figure A.4: Comparison of gufm1 and observatory annual means (SV).
Comparison of the secular variation (SV) predicted by the model gufm1 and annual
means of the secular variation components (X_ = @B(cid:18), Y_ = @B(cid:30) and Z_ = @Br) mea-
(cid:0) @t @t (cid:0) @t
sured at observatories in di(cid:11)erent locations (Fuquene, Mauritius, Niemegk and Tuscon).
Observatory data were obtained from the world data centre. The green triangles are the
observatory annual means of secular variation and the yellow line is the prediction of
gufm1. Only results from the (cid:12)rst observatory site on Mauritius are presented in (b).
257 Appendix A | gufm1
the secular variation observations, but gufm1 again does a good job of (cid:12)tting the trend
in these data.
Comparisonwithdirect,highqualityobservationsmagnetic(cid:12)eldmeasurementsatEarth’s
surfaceshowsthat gufm1is an excellent description of themain (cid:12)eld andits secular vari-
ation, giving con(cid:12)dence in its use in this thesis to analyse patterns of (cid:12)eld evolution.
A.4 Limitations of gufm1
Although gufm1 is our most detailed picture of how the geomagnetic (cid:12)eld has varied in
space and time over the past 400 years, it does have its limitations, in particular those
due to the changing quality and distribution of the available observations with time.
The mapping of surface observations down to the core surface means this variability is
manifested in the changing width of the resolving kernel (the (cid:12)lter through which we
must view the (cid:12)eld at the CMB) with space and time. Backus et al. (1996) discusses
that inferences about the core (cid:12)eld on spatial scales smaller than the area of averaging
kernel peak are impossible, referring to this as the circle of confusion. Jackson (1989)
has produced maps showing how the width of the averaging kernel varies with geograph-
ical position and (cid:12)nds it closely follows the density of observations, so for example in
1882.5 the width of the kernel varies from 11.5 degrees below Europe to 16.5 degrees
in Antarctica where there were no measurements. Animations of gufm1 show the (cid:12)eld
model becomes increasingly complex as time passes, indicating a decrease in the circle of
confusion with time; the smallest features become visible especially in the last 20 years
when satellite data has become available.
Themodelgufm1isultimatelyconstrained(even withperfectdatacoverage) bythelimit
placed on spatial resolution by crustal (cid:12)eld contamination. Unless a sensible method of
removing the crustal (cid:12)eld from the signal can be found, our images of the magnetic (cid:12)eld
at the core surface will always be limited to less than spherical harmonic degree l=12
(around 1800 km at the equator of the core surface), at which point the crustal (cid:12)eld
begins to dominate, and where the (cid:12)eld modelling strategy employed in gufm1 means
that regularisation rather than the (cid:12)t to the data controls the (cid:12)eld structure. In (cid:12)gure
A.5thespatial, sphericalharmonicpowerspectra(as de(cid:12)nedbyLowes (1974)) is plotted
every 50 years to show how the information content of gufm1 changes with time | at
all epochs the power decays monotonically to less than 10(cid:0)3% of the total power by
spherical harmonic degree 10. It is also clear that the more recent epochs contain more
power at lower harmonics than the older epochs. This is because of the requirement
that the models (cid:12)t the data given the error estimates, therefore a larger mis(cid:12)t (hence
smoother model) is permitted when data errors are larger at earlier epochs. The subject
of absolute (cid:12)eld amplitude errors inherent in the core surface (cid:12)eld maps is troublesome.
Gubbinsand Bloxham (1985) producedvariance estimates for their radial magnetic (cid:12)eld
258 Appendix A | gufm1
Figure A.5: Spherical harmonic power spectra from gufm1.
Lowes power spectraof gufm1showingthechange in spectralpower atEarth’s surfaceof
thegeomagnetic(cid:12)eld(cid:12)eldgufm1andillusstrateslimitstoresolutionimposedbydamping
in order to prevent signal contamination by the crustal (cid:12)eld.
models based on the error covariance matrix C and the asymptotic form of the model
e
covariance matrix C (where C(cid:0)1 = (cid:21) S(cid:0)1 +(cid:21) T(cid:0)1 in gufm1). Unfortunately, these
m m S T
errorestimates dependon thechoice of dampingparameters for whichatpresentthereis
noobjective way of determining(theOhmicheating boundsuggests much weaker damp-
ing than is needed to suppress crustal (cid:12)elds | see Backus (1988)). The methodology
employed here solves the existence problem but not the uniqueness problem: as such it
appears di(cid:14)cult to place rigorous bounds on (cid:12)eld amplitude errors.
Despite these limitations, gufm1 is undoubtedly a good model of how the magnetic (cid:12)eld
at the core surface has changed over the past four hundred years. It is consistent with
a huge number of measurements and shows the continuous changes in space and time
expected of a real magnetic (cid:12)eld. Core (cid:13)ow models derived from it can kinematically
account for observed changes in the length of day in the 20th century (Jackson et al.,
1993), an independent con(cid:12)rmation that the models are reasonable.
A.5 Griding of gufm1 for analysis on core surface
In order to carry out a space-time analysis of the magnetic (cid:12)eld features in gufm1, the
spherical harmonic (cid:12)eld model for V was used along with equation (A.8) to evaluate B
r
onaregularlyspacedgridinlatitude, longitudeandtimeonthecoresurface. Thespatial
grid was chosen to be 2 degrees in latitude and longitude and 2 years in time. These are
both much (cid:12)ner scale then the original (cid:12)eld model, ensuring that no information was
lost as a result of the griding procedure.