A Comprehensive Model of The Quiet-Time, Near-Earth Magnetic Field: Phase 3
A Comprehensive Model of The Quiet-Time, Near-Earth Magnetic Field: Phase 3
A Comprehensive Model of The Quiet-Time, Near-Earth Magnetic Field: Phase 3
q
= arccos
_
R
E
+h
q
R
E
+h
A
_
1,2
. (1)
where R
E
=6371.2 km is the mean radius of the Earth and h
q
is al-
titude. Thus,
q
is dened as the latitude of the footprint of an axial
dipole eld line through the apex, and goes to zero at the magnetic
equator at all altitudes. Hence QD coordinates are dened by trac-
ing the eld line upward using DGRF 1980, and downward again
using an axial dipole. This coordinate system is useful in describ-
ing ionospheric and coupling currents that are strongly horizontally
organized by the ambient magnetic eld.
As for the temporal coordinates, there are basicallytworenderings
of time used. The rst is the familiar universal time (UT), which
tracks secular time at Greenwich and is used to describe long-term
and seasonal variations in the elds. The other is magnetic local
time (MLT), t
mlt
, which is dened for an observer as
t
mlt
= (180
+
d.o
d.s
) ,15. (2)
where if the dipole longitude of the observer,
d.o
, and the sub-
solar point,
d.s
, are in degrees, then t
mlt
is in hours. Again, the
tilted dipole from the GSFC(12,83) model is used to assign t
mlt
in
this study. Note that the time angle is now about the dipole z-axis
instead of the geographic z-axis, as in UT, and that a point xed
in MLT is Sun-synchronous. This is useful for describing external
elds whose current systems exhibit diurnal cycles and are also
inuenced by Earth-xed internal elds. For modelling purposes,
diurnal variability is represented by functions expanded in magnetic
universal time (MUT), t
mut
, which is simply the MLT at the dipole
prime meridian, i.e. set
d.o
= 0 in eq. (2).
2 DATA
The accuracy of models derived frominverse problems is intimately
related to the quality of the data being analysed. Undoubtably some
of the best data to date for the purposes of this study come from
the Magsat and POGO satellite missions and from the permanent
magnetic observatories, which will now be discussed.
2.1 Observatories
The CM3 model incorporates two samplings of OHMs fromperma-
nent magnetic observatories. The rst, designated OHM-1AM, are
the OHM values for which the MLTs are closest to 1 am, within 4 h,
chosen fromthe magnetically quietest day of each month, as dened
by the K
p
level, within the interval 19601985. These data were fur-
nished by the World Data Center for Geomagnetismin Copenhagen,
Denmark, and were acquired at a somewhat latter date than the sec-
ond sampling type in order to replace the OAMs rst used in the
model. These offer control of the main eld SV between the Magsat
and POGO mission duration envelopes, which is discussed further
in Section 3.1.
The second, designated OHM-MUL, are chosen from the mag-
netically quietest day of each month during the operational periods
of the POGO (1965 September to 1971 August) and Magsat (1979
November to 1980 April) missions, though data through 1982 are
also included. These data are furnished by the National Geophysical
Data Center (NGDC) in Boulder, CO, with augmented data from
Winch, Faynberg and Singer, and others (priv. comm.). Because the
shortest timescale consideredinthe ionospheric andmagnetospheric
portions of the CM3 model is 6 h (see Sections 3.2 and 3.3), only
OHMs from every other hour are used. Before the OHM-MUL data
set was actually analysed in the CM3 model, it underwent an outlier
rejection phase by visual inspection with respect to a model derived
in the preliminary stages of this study. This process is discussed in
more detail in Section 4.2.
Station breaks, introduced by Langel et al. (1982), are times at
which jumps in the observatory baselines can occur. They are as-
signed by a visual inspection of the time-series of vector eld com-
ponents. These breaks will usually coincide with a physical change
in the nature of the measurements, such as a change in location,
C
2002 RAS, GJI, 151, 3268
36 T. J. Sabaka, N. Olsen and R. A. Langel
equipment or local man-made elds. Thus, each station segment
may be thought to have its own baseline, which is estimated in the
form of a vector bias as described in Section 3.1.
The total measurement counts are listed in the number column
of Table 3 in Section 5, where the OHM-MULs have been divided
into those with dipole colatitude poleward or equatorward of 50
.
The spatial and temporal distributions of the OHM data sets are
shown in Fig. 2. The top panel shows the expected high percentage
of commonality between the two sampling types (triangles). The
bottompanel shows peak station participation during 1965 and 1980
Observatory spatial distribution
180 210 240 270 300 330 0 30 60 90 120 150 180
-90
-60
-30
0
30
60
90
Observatory temporal distribution
0
50
100
150
1960 1965 1970 1975 1980 1985
0
50
100
150
1960 1965 1970 1975 1980 1985
OHM-1AM OHM-MUL
Figure 2. OHM-1AM and OHM-MUL spatial and temporal distributions. The top panel shows observatory locations where only OHM-MUL (squares),
only OHM-1AM (circles), or both (triangles) data sampling types are used in the CM3 model (cylindrical equidistant projection). The bottom panel shows a
histogram of the number of OHM-MUL and OHM-1AM stations contributing data to that particular 1 yr bin across the time span of the CM3 model.
with a slight decrease prior to 1984. Note here that the histograms are
only recording the number of stations, including their breaks, that
provide measurements within a 1 yr bin, and not the total number of
measurements, in which case the OHM-MUL counts would dwarf
those of the OHM-1AM.
2.2 Magsat
The Magsat data sets used in the CM3 model are mainly those
dawn and uncorrected dusk data sets used in deriving the GSFC
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 37
DAWN(6,83) and DUSK(6,83) models of Langel &Estes (1985a),
respectively. Though these data sets are described in detail by Langel
& Baldwin (1991), a brief synopsis of the processing is given here:
the Magsat data were initially screened with the three-hourly K
p
indexbychoosingdata onlyduringperiods when K
p
- 1
following
a periodwhen K
p
2
0
. Vector data polewardof 50
dipole latitude
were excluded to minimize the effects of FACs and ionospheric
currents in the auroral regions. Scalar intensity data were retained at
all latitudes. A data selection algorithm was then applied separately
for both the dawn and dusk data within a 20 nT D
st
level for the
time intervals of 1979 November and December; 1980 January and
February, and 1980 March and April. The objective of the algorithm
was to obtain a uniform data distribution in both time and space.
Finally, after elimination of outliers with respect to the GSFC(9,80)
model of Langel et al. (1982), passes from slightly more disturbed
times were added to sparse regions in order to improve geographic
coverage. These data sets will be referred to simply as the Magsat
dawn and Magsat dusk data sets.
The CM3 model represents a signicant advancement in how
elds of external origin, such as those from ionospheric coupling
and EEJ currents, are parametrized. It is then desirable to include
high (poleward of 50
, correspondingtothose already
providing scalar measurements in the CM3 Magsat dawn and dusk
data sets. These new data sets will be referred to as the Magsat
polar dawn and Magsat polar dusk data sets.
Before the Magsat dawn, polar dawn, dusk and polar dusk data
sets were actually analysed in the CM3 model, they underwent an
additional episode of renement via the rejection of outliers with
respect to a model derived in the preliminary stages of this study.
This process is discussed in more detail in Section 4.2, and the
resulting measurement counts for these data sets are listed in the
number column of Table 3 (see Section 5).
2.3 POGO
The bulk of the parent POGO data sets from which those used
in the CM3 model were extracted were also used to derive the
POGO(2,72) eldmodel (Langel et al. 1980), withadditional OGO-
6 data from 1969 to 1971 for magnetically quiet to moderately quiet
times. These parent data sets, described in detail by Langel & Bald-
win (1991), were found to have an uneven distribution in local time
and to be somewhat larger than necessary for this study. Therefore,
these data sets were decimated in order to achieve a more uniform
geographic and magnetic local time distribution, and a more man-
ageable size. Furthermore, the same outlier rejection phase was ap-
plied to POGO data as to Magsat data, details of which are found in
Section 4.2. The resulting data sets will be referred to collectively
as the POGO decimated data set. The MLT distribution of this
data set is shown in the bottom panel of Fig. 3 as a histogram of
the number of measurement positions falling within 1 h MLT bins.
Of the total 22 685 positions, most (1509) fall between 1:00 and
2:00 am and the least (390), unfortunately, fall between 12:00 and
1:00 pm.
Historically, the mechanism used for the POGO decimated data
set did not admit entire satellite tracks. Much of the external eld
current systemis transient and, while this distribution probably gives
a broad sampling of those variations, it may not give coherent data
tracks across these patterns. Because of this, it is felt that some
advantage might be gained by incorporating some individual passes
of data. Accordingly, a selection of data from typical passes from
quiet periods has been made. These data were also put through the
outlier rejection phase, which resulted in 6754 measurements from
170 passes. This data set is referred to as the POGO pass data set.
The angular positions of the pass loci are shown in the top panel
of Fig. 3, and a histogram of the number of passes that cross the
geographic equator within 1 h MLT bins is shown in the middle
panel. The spatial and temporal distributions appear to be sufcient
to sample most of the Sq and EEJ features of the ionospheric current
systems.
3 PARAMETRI ZATI ON
OF FI ELD S OURCES
As alluded to earlier, the ultimate worth of the CM3 model lies in
its ability to properly describe as much of the intended signal as
possible in the data. These data, however, are limited in their spatial
and temporal sampling of the eld and are contaminated at some
level by systematic and random error processes. An efcient model
parametrization will take these limitations into account in order to
achieve optimal results. The parametrizations used in this study are
now described by source.
3.1 Core and lithospheric elds
The current systems responsible for both the core and lithospheric
magnetic elds lie entirely belowthe regions sampled by permanent
observatories andsatellites. Therefore, these elds maybe expressed
as gradients of internal potential functions of the form
V
cl
(t. r) =
_
a
Nmax
n=1
n
m=0
_
a
r
_
n+1
m
n
(t )Y
m
n
(. )
_
(3)
with
Y
m
n
(. ) = P
m
n
(cos ) exp i m. (4)
where a is the mean radius of the Earth (6371.2 km), r is the position
vector, and Y
m
n
and P
m
n
are the Schmidt normalized surface spherical
harmonic and associated Legendre function of degree n and order
m, respectively (e.g. Langel 1987). The {} operation takes the
real part of the expression only. Hence,
m
n
are unique complex
expansion coefcients, also known as Gauss coefcients. They are
related to the usual real Gauss coefcients g
m
n
and h
m
n
according to
m
n
= g
m
n
i h
m
n
. Typically, terms in eq. (3) have been retained only
up to a degree truncation level, N
max
, that is justied by the data, or in
the case of satellites, up to the degree at which it is believed that the
lithospheric eldbegins to dominate the series (the main eld), taken
by Langel & Estes (1982) to be 13. Spherical harmonic models of
lithospheric elds have been derived fromdata with estimates of the
main, magnetospheric and ionospheric elds removed. Such models
indicate that noise becomes dominant somewhere between N
max
=
60 and 70 (Ravat et al. 1995). That noise, however, also reects an
inaccurate estimation of the other elds. One of the intentions of
the present study is to examine whether a combined model can be
more effective in accurately representing the lithospheric eld, and
so N
max
= 65 in this study.
3.1.1 Main eld secular variation
The main eld is mostly dominated by contributions that originate
in the core dynamo region of the Earth. The consequence of this
action is that the main eld varies with timescales of the order of
C
2002 RAS, GJI, 151, 3268
38 T. J. Sabaka, N. Olsen and R. A. Langel
POGO pass spatial distribution
180 210 240 270 300 330 0 30 60 90 120 150 180
-90
-60
-30
0
30
60
90
POGO pass equatorial crossing MLT distribution
0
4
8
12
0 4 8 12 16 20 24
POGO decimated MLT distribution
0
500
1000
1500
0 4 8 12 16 20 24
Figure 3. POGO satellite spatial and temporal distributions. The top panel shows the angular positions of the pass loci for the POGO pass data set (cylindrical
equidistant projection). The middle panel shows a histogram of the number of passes from the POGO pass data set that cross the geographic equator within
a particular 1 h MLT bin. The bottom panel shows a histogram of the number of measurement positions in the POGO decimated data set that fall within a
particular 1 h MLT bin.
centuries. Variations with timescales shorter than a year or two are
effectively screened by nite mantle conductivity, and as a result, are
not detectable at the Earths surface. Sabaka et al. (1997) modelled
the longer timescale variations of the main eld by parametrizing its
SVwith cubic B-splines, which are described in detail by Schumaker
(1981).
This approach is adopted here using cubic B-splines and the 1980
Magsat epoch for all
m
n
(t ). Because both POGO and Magsat satel-
lite data will be included in the analysis, the time span of the model
is chosen to be 19601985, allowing for an extension at either end of
the missions. Given that the spline knot set must be dened over the
entire time envelope, and choosing an equispaced knot distribution
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 39
at 2.5 yr intervals for all
m
n
(t ), this results in 13 cubic B-splines per
m
n
(t ). Including the offset,
m
n0
, this gives a total of 14 parameters
describing the temporal behaviour of each
m
n
(t ).
Assuming that the temporal variation of terms with n >13 is
negligible for the present model, the nal working expression for
the core and lithospheric potential is given by
V
cl
(t. r) =
_
a
13
n=1
n
m=0
13
q=0
_
a
r
_
n+1
m
nq
Y
m
nq
(t. . )
+ a
65
n=14
n
m=0
_
a
r
_
n+1
m
n
Y
m
n
(. )
_
. (5)
with
Y
m
nq
(t. . ) =
_
_
Y
m
n
(. ) for q = 0.
Y
m
n
(. )
_
t
1980
b
q
() d for q > 0.
(6)
where b
q
() is the qth cubic B-spline of the expansion. The number
of real coefcients in the main and lithospheric eld model expan-
sions are 2730 and 4160, respectively, which gives a total of 6890
real coefcients in this portion of the model.
3.1.2 Observatory biases
Permanent magnetic observatories are extremely sensitive to the
short-wavelength eld of the lithosphere; even at N
max
=65 there is
undoubtedly deviation owing to an unmodelled high-degree litho-
spheric signal. Hence, this signal is point-sampled by the observato-
ries and can be represented by a local vector bias, B
bias
, as introduced
by Langel et al. (1982). There are 427 individual observatory time-
series in the combined OHM-1AM and OHM-MUL data sets for
which a full vector bias is estimated. Owing to data quality issues,
only the horizontal data are used and bias components estimated
for the Amberley II OHM-MULs and the San Fernando OHM-
1AMs and OHM-MULs. Similarly, only the vertical data are used at
Amatsia, Gornotayezhnaya, Norway Station and Sheshan, and
hence, only this component of the biases is estimated. Note that
although data from a particular station may be present in both the
OHM-1AM and OHM-MUL data sets, separated biases are esti-
mated. This has little effect, however, and was done mainly to ex-
pedite the replacement of the OAMs by the OHM-1AMs. It will
be resolved in future versions of the model. The biases account for
1296 real coefcients in the model.
3.2 Ionospheric eld
The morphology of geomagnetic variations produced by the iono-
spheric dynamo is relatively xed in MLT. However, within this
basic morphology there is considerable variation depending upon
other phenomena such as the season, the solar cycle, the inuence of
the Earths main eld, etc. The ionospheric Sq currents are conned
to the E-region dynamo layer (100140 kmaltitude), with maximum
current density at about 110 km altitude. Since the thickness of this
layer (-50 km) is small compared with the distance to the point
of observation (>100 km), the true ionospheric current can be well
approximated by an equivalent sheet current at h =110 km altitude,
which produces the same magnetic eld outside the current distri-
bution. Because these elds vary with time, there are corresponding
induced currents in the Earth, with attendant elds. Since the mea-
surements used are not acquired in the regions where the source
currents ow, the elds can be represented by gradients of potential
functions. The basis functions representing the ionospheric and as-
sociated induced potentials are best understood as being built from
a set of elemental potential functions reecting a single spatial
harmonic modulated by single seasonal and diurnal periods. For the
region a r a +h these have the form
V
m
nsp
(t. t
mut
. r) =
_
m
nsp
S
m
nsp.i
(t. t
mut
. r) +c
m
nsp
S
m
nsp.e
(t. t
mut
. r)
_
.
(7)
and for the region r > a +h these have the form
(V
)
m
nsp
(t. t
mut
. r) =
__
m
nsp
+(c
)
m
nsp
_
S
m
nsp.i
(t. t
mut
. r)
_
. (8)
with
S
m
nsp.e
(t. t
mut
. r) =a
_
r
a
_
n
P
m
n
(cos
d
)
exp i (m
d
+
s
st +
p
pt
mut
). (9)
S
m
nsp.i
(t. t
mut
. r) =a
_
a
r
_
n+1
P
m
n
(cos
d
)
exp i (m
d
+
s
st +
p
pt
mut
). (10)
where the fundamental seasonal angular frequency is
s
= 2 rad
yr
1
with associated wavenumber s and the fundamental diurnal
angular frequency is
p
= 2,24 rad h
1
with associated wave-
number p. Hence, c
m
nsp
, (c
)
m
nsp
and
m
nsp
are unique complex expan-
sion coefcients of the external and internal ionospheric and the
internal induced potentials, respectively, having a single spatial har-
monic as prescribed by n and m, which oscillates on two timescales
as prescribed by s and p, and propagates in one direction as pre-
scribed by the relative signs of s, p and m. However, the ionospheric
potentials above and below the sheet source are by no means inde-
pendent because the radial components of the resulting elds are
continuous across the sheet. This results in the linear relationship
(c
)
m
nsp
=
_
n
n +1
__
a +h
a
_
2n+1
c
m
nsp
. (11)
such that eq. (8) may be rewritten in terms of c
m
nsp
.
3.2.1 Induction
At this point one could dene an ensemble of V
m
nsp
as the working
formof the ionospheric and associated induced potentials. However,
incorporating a few conditions gained from some basic physical
insights can drastically reduce the number of free parameters in this
part of the model. The rst insight that can be made is that
m
nsp
is, in
general, not independent of c
m
nsp
. The nature of this dependence is
related to the conductivity structure of the crust and upper mantle,
which leads to the linear relationship
l
ksp
=
n.m
q
lm
knsp
c
m
nsp
. (12)
where q
lm
knsp
is an element of the complex matrix representation, Q,
of the transfer function between the driving ionospheric signal and
the driven induced signal at similar frequencies (Schmucker 1985;
Olsen 1999). For a general 3-D mantle conductivity, the Q matrix
is dense. However, for purposes of this model the 1-D conductivity
distribution, i.e. depending only on radius, of Olsen (1998) has been
adopted. It is a four-layer conductivity model derived from Sq and
D
st
data at selected European observatories. Consequently, Qis now
diagonal and its elements depend only upon n and the frequency f.
This means that one external coefcient induces only one internal
coefcient, i.e.
m
n
( f ) = q
mm
nn
( f )c
m
n
( f ) (where the f dependence
C
2002 RAS, GJI, 151, 3268
40 T. J. Sabaka, N. Olsen and R. A. Langel
now replaces the s and p indices). The mechanics of calculating Q
from a given conductivity model is beyond the scope of this paper,
but the interested reader is referred to Parkinson & Hutton (1989).
The value of |q
mm
nn
( f )| typically decreases with decreasing fre-
quency, and therefore elds induced from higher frequencies gen-
erally dominate those from lower ones. Therefore, for simplicity, a
single frequency (the highest of the diurnalseasonal modulation)
may be assigned to the matrix elements operating on c
m
nsp
, which
leads to the assignment rules
q
mm
nn
( f ) =
_
_
0 for p = 0 and s = 0.
q
mm
nn
(0) for p = 0 and |s| > 0.
q
mm
nn
( p) for p > 0.
(13)
Note that f =0 is formally used to designate very long-period
rather thananabsence of oscillation, inwhichcase q
mm
nn
( f 0) =0.
This is indeed true for the purely seasonal versus diurnal oscillations
included in the model. For deriving the very long-period response
q
mm
nn
( f 0), the mantle is assumed to be an insulator in the region
a r a and superconducting in r - a . For this study,
a value of =1000 km is used, corresponding to induction with
periods longer than 1 week or so.
3.2.2 Quasi-dipole symmetry
A second insight is that many ionospheric phenomena are naturally
organized with respect to the geometry of the Earths magnetic eld
owing to its inuence on the motion of charged particles. Conse-
quently, ionospheric conductivity is highly anisotropic, resulting in
values that are so high parallel to the eld that the magnetic eld
lines are nearly equipotential lines. Therefore, it is often convenient
to work in a coordinate system that is aligned with the magnetic
eld, such as the QD coordinates proposed by Richmond (1995). If
one could exploit the symmetries of such a coordinate system, then
great savings would be realized in the number of free parameters
required. For example, modelling the EEJ with spherical harmonics
expanded in dipole coordinates requires both high-degree and high-
order terms because of the undulation of the EEJ in
d
. However, the
EEJ is always located at
q
= 90
n=1
min (n.Mmax)
m=min (n.Mmax)
_
d
lm
kn
_
P
m
n
(cos
d
) exp i m
d
. (14)
where
q
=
q
(r.
d
.
d
),
q
=
q
(r.
d
.
d
), d
lm
kn
= (d
lm
kn
)
, and
the asterisk denotes complex conjugation. The regression coef-
cients, d
lm
kn
, are determined by a standard spherical transform,
where N
max
and M
max
are chosen such that sufcient convergence
is achieved. The required temporal dependence is achieved by
simply multiplying by the appropriate complex exponential, e.g.
exp i (
s
st +
p
pt
mut
).
If one settles for the radial dependences found in eqs (9) and
(10), then the new functions will satisfy Laplaces equation in
dipole coordinates for external and internal sources, respectively,
and will exhibit the desired QD symmetry on a given sphere. Since
h =r a =110 km is the approximate height of the ionospheric
E-region current system, then this should be a satisfactory choice.
Hence, the potentials may be expressed in terms of
T
l
ksp.e
(t. t
mut
. r) =
n.m
_
d
lm
kn
_
_
a
a +h
_
n1
S
m
nsp.e
(t. t
mut
. r). (15)
T
l
ksp.i
(t. t
mut
. r) =
n.m
_
d
lm
kn
_
_
a +h
a
_
n+2
S
m
nsp.i
(t. t
mut
. r). (16)
where the summations over n and m are understood from eq. (14).
Clearly, it follows that
T
l
ksp.e
r=a+h
=T
l
ksp.i
r=a+h
=(a +h)P
l
k
(cos
q
) exp i (l
q
+
s
st +
p
pt
mut
).
(17)
Aglobal plot of {T
3
45.0.0.e
}|
r=a+h
is shown in Fig. 4 which illustrates
the QD geometry on the sphere where the ionospheric currents are
assumed to ow.
The previous development suggests a new set of elemental po-
tential functions, analogous to those of eqs (7) and (8), the basis of
which consists of T
l
ksp.e
and T
l
ksp.i
, which indicates that imposing the
QD symmetry at r = a +h imposes linear constraints on the origi-
nal expansion coefcients. Thus, by combining eqs (11), (12), (15)
and (16), one arrives at expressions for the expansion coefcients
of the original ionospheric and induced potentials, given by
c
m
nsp
=
k.l
d
lm
kn.e
c
l
ksp
. (18)
(c
)
m
nsp
=
k.l
z
lm
kn
c
l
ksp
. (19)
m
nsp
=
k.l
f
lm
knsp
c
l
ksp
. (20)
where
d
lm
kn.e
=
_
a
a +h
_
n1
d
lm
kn
. (21)
z
lm
kn
=
_
n
n +1
__
a +h
a
_
n+2
d
lm
kn
. (22)
f
lm
knsp
= q
mm
nnsp
d
lm
kn.e
. (23)
in terms of the QD ionospheric expansion coefcients, c
l
ksp
.
3.2.3 Solar activity
Ionospheric contributions also depend upon solar activity. In a paper
on the variability of geomagnetic daily variations with solar activity,
Olsen (1993) estimates the proportionality between the coefcients
of a spherical harmonic analysis of the variations and sunspot num-
ber. Solar ux, however, is probably a better parameter for describing
the short-term variability of solar activity of the kind encountered
in this study. Therefore, the QD ionospheric expansion coefcients
are redened with a dependence upon the solar radiation ux index,
F
10.7
, as
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 41
Figure 4. The function {T
3
45.0.0.e
} evaluated on the sphere r = 6481.2 km (cylindrical equidistant projection). The reference model used to dene the QD
coordinate system is the DGRF 1980 model (IAGA Division I Working Group 1, 1981). The
q
and
q
coordinate lines are also shown in 30
increments.
c
l
ksp
= c
l
ksp
(1 + NF
10.7
) . (24)
The proportionality N is not solved for, but rather an indepen-
dently estimated a priori value is used, which is assumed to be
equal for all coefcients. Thus, N was estimated using the tech-
nique of Olsen (1993), resulting in a value of N =14.85
10
3
(10
22
Wm
2
Hz
1
)
1
. This means that increasing solar ux
inates the whole ionospheric current system (and induced con-
tributions) without changing its shape. In this study, annual av-
erages of absolute F
10.7
values were used for all data sets ex-
cept the OHM-1AMs, which were inadvertently assigned monthly
averages of adjusted F
10.7
values. Adjusted values are corrected
for the changing SunEarth distance, while absolute values have
an additional correction for uncertainties in antenna gain and for
waves reected from the ground, which amounts to multiplying
the adjusted value by 0.9. For more information on the F
10.7
index, see http://web.ngdc.noaa.gov/stp/SOLAR/FLUX/ux.html.
However, this is felt to be a minor problem since nightside iono-
spheric E-region conductivity is so low; a property that is enforced
in this model (see Section 4.1.2).
3.2.4 Basis function selection
Given the temporal distribution of the magnetic measurements, the
rst four diurnal harmonics can probably be resolved. This corre-
sponds to p =1. . . . . 4, which is to say, the 24, 12, 8 and 6 h periods.
The p =0level will alsobe includedfor reasons tobe discussedlater.
Instead of using all QD longitudinal wavenumbers l, it is suitable to
focus on local time modes ( p = l) plus slightly faster and slightly
slower travelling modes. This limits l to a tight bandwidth, L, centred
on p. In this study, L =1. To obtain similar QDlatitudinal resolution
across diurnal periods, the maximum k should be a constant offset,
K, from l. It is also desirable to have different resolution levels for
local ( p = l) versus non-local ( p = l) time modes, in which case
K becomes a function of p l (i.e. K = K( p l)). In this study,
they happen to be the same, with K(0) = K(= 0) = 40, chosen in
hopes of resolving the EEJ. Note that this general selection scheme
for p, k and l has been used in previous studies to produce global
eld maps from Sq currents (e.g. Malin & Gupta 1977; Schmucker
1999).
Seasonal variation is, to rst order, a function of the angle be-
tween the EarthSun line and the Earths rotation axis. However,
it is also inuenced by the Earths magnetic eld, for which the
dipole portion alone wobbles about the rotation axis daily. This, to-
gether with comingling of other effects, such as the solar radiation
ux with its own associated periodicities, makes for a very complex
picture indeed. It is believed, however, that the annual and semi-
annual periods will still dominate the seasonal spectrum, and so
these are the modes that will be considered in this model. Further-
more, both eastward and westward counterparts of these modes will
be included. All this translates into a range of s =2. . . . . 2 for the
model.
The number of columns of the matrices associated with ele-
ments d
lm
kn.e
, f
lm
knsp
and z
lm
kn
are specied by the s, p, k and l ranges.
The number of rows is determined by N
max
and M
max
in eq. (14),
which are taken to be 60 and 12, respectively, based upon values
of max (k) =45 and max (|l|) =5. This gives a total of 1368 real
regression coefcients per T
l
ksp.e
.
C
2002 RAS, GJI, 151, 3268
42 T. J. Sabaka, N. Olsen and R. A. Langel
Thus, the pair of expressions for the ionospheric and associated
induced potentials for the regions a r a + h and r > a + h
may be written as
V
ion
(t. t
mut
. r) =
_
2
s=2
4
p=0
p+1
l=p1
|l|+40
k=max (1.|l|)
c
l
ksp
60
n=1
min (n.12)
m=min (n.12)
__
d
lm
kn.e
_
S
m
nsp.e
(t. t
mut
. r)
+
_
f
lm
knsp
_
S
m
nsp.i
(t. t
mut
. r)
_
_
_
_
. (25)
V
ion
(t. t
mut
. r) =
_
2
s=2
4
p=0
p+1
l=p1
|l|+40
k=max (1.|l|)
c
l
ksp
60
n=1
min (n.12)
m=min (n.12)
__
z
lm
kn
_
+
_
f
lm
knsp
_
_
S
m
nsp.i
(t. t
mut
. r)
_
_
_
. (26)
The quadruple summation over s, p, l and k implies 5520 real coef-
cients. This is about six times less than the number of parameters
that would be needed if QD symmetry were not considered.
In summary, Sq currents are organized spatially and temporally by
the interaction of the ambient magnetic eld and the Sun. Hence, to
model the associated elds efciently this geometry is exploited in
the design and selection of the basis functions. This reduces to a ju-
dicious choice of linear combinations of the elemental functions so
as to basically followQDsymmetry in space and Sun-synchronicity
in time, allowing for some non-local time interaction and solar ac-
tivity dependence. Proper coupling of primary and induced elds
and of elds above and below the equivalent sheet current gives a
representation of these elds in terms of a single set of expansion
coefcients, c
l
ksp
.
3.3 Magnetospheric eld
The eld of the magnetosphere is dominated by features that vary
with ring-current intensity, season and solar wind parameters, but
also contains features that are relatively xed in MLT. These time-
varying elds also induce currents and resultant secondary elds in
the conductive portions of the crust and mantle. Because the major
contributingsources of the magnetospheric eldare the magnetotail,
magnetopause and ring-current complexes, which lie well outside
the sampling region, the eld may be represented by the gradient
of a potential function. In fact, the form and development of this
function will closely parallel that of the ionosphere and associated
induced elds for the region a r a + h. Therefore, one may
anticipate a nal working form for the magnetospheric potential
very similar to eq. (25), and expedite its development by discussing
only the deviations from the ionospheric case.
In addition to the regular daily and seasonal periodicities, a mod-
ulation of the magnetospheric elds with the D
st
index, given by
Langel & Estes (1985a) as
c
m
1sp
= j
m
1sp
+j
m
1sp.Dst
D
st
(t ). (27)
is used here, where the D
st
index is in units of nT and is tabulated at
hourly intervals as a function of UT (see http://swdcdb.kugi.kyoto-
u.ac.jp/). Note that this relationship is adopted only for the dipole
terms (n =1), andthat the temporal variabilityof D
st
(t ) is modulated
by both seasonal and diurnal oscillations to help describe any local
time asymmetries.
3.3.1 Induction
Induced elds may be treated just as they were in the ionospheric
eld development section. The major difference is the inclusion of
basis functions that are nowdependent upon the seasonally and diur-
nally modulated D
st
(t ) index. Hence, for magnetospheric induction,
eq. (13) becomes
q
mm
nn
( f ) =
_
_
0 for p = 0. s = 0. and noD
st
(t ).
q
mm
nn
(0) for p = 0. s = 0. and D
st
(t ).
q
mm
nn
(0) for p = 0 and |s| > 0.
q
mm
nn
( p) for p > 0.
(28)
Magnetospheric contributions that varyonlywith D
st
(t ) are assumed
to contain mostly signals with a period of a few days, at least during
magnetically quiet days, hence, the use of a q
mm
nn
(0) based upon
= 1000 km.
3.3.2 Basis function selection
At the source region for the magnetospheric current systems the
Earths magnetic eld is more dipole-like compared with iono-
spheric current systems, and therefore, there is no benet in using
QD coordinates to characterize magnetospheric sources. Again, l
is chosen to reside in a narrow band about p, as specied by L,
and maximum k is at a constant offset, K, from l to preserve lat-
itudinal resolution levels across p. K is a function of local versus
non-local time, as before. A seasonal variation of s =2. . . . . 2, a
diurnal variation of p =0. . . . . 5, and L =1 are also used. In con-
trast with the ionosphere, the latitudinal resolution level, which is
much less for the magnetosphere and different for local and non-
local time modes, is described by K(0) =5 and K(= 0) =3. It is
expected that a signicant number of coefcients will be negligi-
ble at these truncation levels. In an investigation of geomagnetic
daily variations as predicted by the Tsyganenko model of the mag-
netosphere, Olsen (1996) concluded that the only non-negligible
coefcients are found within the limits of p = 0. . . . . 2, L = 1,
and K(0) =3 and K(=0) =1. However, since this is the rst time
ionospheric and magnetospheric parametrizations of this type have
been coestimated, these limits are chosen more liberally.
Thus, the magnetospheric and associated induced potentials may
be written as
V
mag
(t. t
mut
. r) =
_
2
s=2
5
p=0
p+1
l=p1
|l|+K( pl)
k=max (1.|l|)
j
l
ksp
_
S
l
ksp.e
(t. t
mut
. r)
+
_
q
ll
kksp
_
S
l
ksp.i
(t. t
mut
. r)
_
+
2
s=2
5
p=0
p+1
l=p1
1
k=max (1.|l|)
j
l
ksp.Dst
D
st
_
S
l
ksp.e
(t. t
mut
. r)
+
_
q
ll
kksp.Dst
_
S
l
ksp.i
(t. t
mut
. r)
_
_
_
_
. (29)
Note the s, p and D
st
indexing on the Qelements, which is consistent
with eq. (28). The summations entail 800 real coefcients.
3.4 Fields from ionospheric coupling currents
In reality, ionospheric currents do not ow in isolated shells, but are
coupled to the magnetospheric and ionospheric E-region currents
at the geomagnetic conjugate point through coupling currents that
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 43
ow along the eld lines of the Earths magnetic eld. This means
that the Magsat samplingregionbetweena+350anda+550kmwill
be penetrated by F-region current for which the associated magnetic
eld, B, will not be curl-free, and hence, will not be expressible as
the gradient of a potential. It should also be noted that in general
only vector samples, as opposed to scalar ones, can detect these
elds, since B is almost always perpendicular to the main eld.
Hence, only vector measurements from the Magsat satellite will be
considered here.
Consider the toroidalpoloidal decomposition of B. Backus
(1986) showed that for a shell with thickness that is small in com-
parison with its mean radius, the poloidal part goes to zero. This
is the case for the Magsat sampling shell, and so B is considered
purely toroidal such that
B(t. r) = r+(t. r). (30)
with associated current density
j
0
J(t. r) = r+(t. r). (31)
where + is the toroidal scalar function or stream function. In order
for B to be unique, the mean value of + over the sphere must be
zero. If + is expanded in surface spherical harmonics as
+(t. r) =
_
n.m
m
n
(t. r)P
m
n
(cos
d
) exp i m
d
_
. (32)
then this requires
0
0
(t. r) = 0.
As for the radial dependence of
m
n
(t. r), Olsen (1997a) consid-
ered a 1,r dependence on all coefcients, which leads to a purely
radial poloidal current density in eq. (31). When applying this to
Magsat data, the assumption of radial currents is reasonable for the
FACs in polar latitudes as well as for a very narrow band at the
dip equator where the meridional current system ows in the ra-
dial direction. However, the assumption fails at middle latitudes.
Nevertheless, this assumption is also adopted here, which leads to
+(t. r) =
_
n.m
m
n
(t )S
m
n.j
(r)
_
. (33)
where
S
m
n.j
(r) =
_
a
r
_
P
m
n
(cos
d
) exp i m
d
. (34)
3.4.1 Quasi-dipole symmetry
The F-region conductivity structure is also highly aligned with the
magnetic eld, which suggests the use of QD symmetry. For exam-
ple, the meridional coupling currents of the EEJ showa strong radial
upwelling along
q
= 90
n.m
_
d
lm
kn
_
_
a +h
a
_
S
m
n.j
(r). (35)
It follows that
T
l
k.j
r=a+h
= P
l
k
(cos
q
) exp il
q
. (36)
The d
lm
kn
regression coefcients would be slightly different from
those introduced earlier since these now reect the QD symmetry at
r = a+450 km. As before, the newbasis imposes linear constraints
on the original expansion coefcients such that
m
n
(t ) =
k.l
d
lm
kn.j
l
k
(t ). (37)
where
d
lm
kn.j
=
_
a +h
a
_
d
lm
kn
(38)
and
l
k
(t ) are the new QD coefcients.
3.4.2 Time dependence
As in the E-region, the temporal variation of the magnetic elds from
the coupling currents will have strong local time modes, which are
modulated by both interactions with the main eld and by signi-
cant seasonal effects. As before, one could introduce the variation
through factors of the form exp i (
s
st +
p
pt
mut
). However, only
Magsat data are sensitive to these parameters, and this data has
some severe limitations. First, the mission duration was approxi-
mately from 1979 November to 1980 April, about six months. If
one is interested in both annual and semi-annual seasonal periods,
then it is unlikely that both a phase and amplitude can be resolved.
Assuming, however, that the maximum annual amplitudes occur at
the solstices, then the phase may be xed, leaving only the am-
plitude to be determined. If it is also assumed that maximum semi-
annual amplitudes occur at the solstices andequinoxes, thenits phase
may also be xed. This leads to admissible basis functions of the
form
T
l
ks.j
(t. r) =
n.m
_
d
lm
kn.j
_
S
m
ns.j
(t. r). (39)
where
S
m
ns.j
(t. r) = S
m
n.j
(r) cos
s
st . (40)
with associated coefcients
l
ks
.
Secondly, the Magsat data are sampled at only two local times,
dawn and dusk. The Nyquist frequency for this sampling rate would
be the diurnal ( p =1) frequency; all higher-frequency harmonics
would be aliased. Therefore, a continuous local time analysis is
prohibited at the periods of interest. Consider, however, that over a
period of several days to weeks a high-density distribution of both
Magsat dawn and dusk passes may be obtained over all longitudes.
This suggests that one model the behaviour of particular local times
separately as a function of geographic position and season.
3.4.3 Basis function selection
Separate dawn and dusk contributions are used, each with a sea-
sonal wavenumber range of s = 1. 2. The k and l ranges for these
contributions are less complicated to specify than for the E-region,
since they reect the spatial resolution level in an Earth-xed mode.
Consequently, the latitudinal and longitudinal resolution levels are
independent of one another. It is anticipated that most features of in-
terest will lie in relatively thin, elongated QD latitude bands, which
implies a high K
max
and a low L
max
. If the eld from the EEJ cou-
pling currents has a half-wavelength of about 5
at
q
=90
, then
K
max
= 40 should sufce. For this study, L
max
=4. Given these val-
ues of K
max
and L
max
, an expansion with N
max
= 60 and M
max
=12
is considered sufcient, i.e. a total of 1368 real regression coef-
cients per T
l
ks.j
. Thus, the stream function may be expressed as
+(t. r) =
_
2
s=0
40
k=1
min (k.4)
l=0
l
ks
60
n=1
min (n.12)
m=min (n.12)
_
d
lm
kn.j
_
S
m
ns.j
(t. r)
_
. (41)
C
2002 RAS, GJI, 151, 3268
44 T. J. Sabaka, N. Olsen and R. A. Langel
Table 1. Number of parameters describing each eld
source in the CM3 model.
Field source Number of parameters
Observatory biases 1296
Core/lithosphere 6890
Ionosphere 5520
Magnetosphere 800
Coupling currents 2088
Total 16 594
The summations entail 1044 real coefcients, so 2088 for both dawn
and dusk. This is about four times fewer than the number of param-
eters if QD symmetry were not considered.
3.5 Observed elds
The expressions for predicting the magnetic eld as seen by Magsat,
POGO and OHMs are given by
B
Magsat
=
g
V
cl
R
gd
[
d
(V
ion
+ V
mag
)
d
r+]. (42)
B
Magsat
=
g
V
cl
R
gd
d
(V
ion
+ V
mag
). (43)
B
POGO
=
g
V
cl
R
gd
d
(V
ion
+ V
mag
). (44)
B
OHM
= R
es
[
g
V
cl
+R
gd
d
(V
ion
+ V
mag
)] +B
bias
. (45)
where V
cl
, V
ion
, V
ion
and V
mag
are the potential functions for the elds
of the core and lithosphere, the ionospheric E-region for a r
a +h and r > a +h, and the magnetosphere, respectively.
g
and
d
are the gradient operators in geographic and dipole spherical
polar coordinates, respectively. R
gd
is a rotation matrix from the
local dipole spherical to the local geographic spherical basis. R
es
is a rotation matrix from the local geographic spherical to the local
geographic ellipsoidal basis, where the IAU ellipsoid (International
Astronomical Union 1966) is used (Barton 1997). A summary of
the number of parameters describing each eld source is given in
Table 1.
4 ES TI MATI ON OF MODEL
PARAMETERS
Having essentially dened the forward problem for this study in
eqs (42) (45), the model parameters were found by solving the in-
verse problemvia the weighted least-squares method by minimizing
the cost function
L(x) = [d a(x)]
T
W(d a(x)) +
k
i =1
i
(x x
a
i
)
T
A
i
(x x
a
i
)
(46)
with respect to the model parameter vector x, where d is the mea-
surement vector, a(x) is the model prediction vector (eqs 4245),
x
a
i
is the preferred parameter vector from the ith constraint system,
W and A
i
are the uncalibrated inverse covariance matrices of the
data noise vector and the deviation vector from x
a
i
, respectively,
i
is the damping parameter fromthe ith constraint system, and k is the
number of constraint systems. The Gauss iterative method was used
to estimate the optimal model. Its nth step is given by (Tarantola &
Valette 1982)
x
n+1
=x
n
+
_
A
T
n
WA
n
+
k
i =1
i
A
i
_
1
_
A
T
n
W(d a(x
n
)) +
k
i =1
i
A
i
(x
a
i
x
n
)
_
. (47)
where x
n+1
is the estimated model parameter vector after n steps
and A
n
is the Jacobian of a(x) at x
n
.
The uncalibrated error-covariance matrix of the nal parameter
estimate, x, is
C
x
=
_
A
T
WA +
k
i =1
i
A
i
_
1
. (48)
where A is the Jacobian at x. An unbiased estimate, s
2
, of the data
mist is given by (Toutenburg 1982)
s
2
=
[d a( x)]
T
W(d a( x))
N tr[R]
. (49)
where
R = C
x
A
T
WA. (50)
and where tr[] is the trace operator and N is the number of obser-
vations. The resolution matrix R acts as a lter through which the
true model state is seen and the trace of which gives the number of
model parameters resolved by the data (Langel 1987). Thus, s
2
is a
measure of how well the model ts the weighted data per degree of
freedom(DOF), and should be approximately unity if the weighting
is correct (Bloxham et al. 1989). One may now obtain an unbiased
estimate of the calibrated parameter error-covariance matrix by sim-
ply multiplying C
x
by s
2
. In this sense a calibrated error-covariance
matrix reects how well the observations can actually be t.
4.1 Regularization and a priori information
All of the constraints considered in this study are quadratic norms
in x, where A
i
is the associated norm matrix, and these will now be
discussed.
4.1.1 Smooth main eld secular variation
In order for CM3 to remain consistent throughout its temporal do-
main the main eld SV must behave in a plausible manner. This is
a critical issue because of the restricted distribution of the data sets:
Magsat and POGOprovide small-scale spatial information, but only
for short durations, while OHMs cover the entire domain, but only
provide broad-scale spatial information. Thus, for CM3 the class
of admissible SV models is restricted by imposing temporal and
spatial smoothing (least-structure) constraints. Since the notion of
smoothing usually implies a minimization of curvature, an attempt
was made to smooth the second-order time and space derivatives
of the radial eld SV,
B
r
, through time over the CMB. The spatial
derivatives actually take the form of the surface Laplacian. How-
ever, it was found that the second-order time derivative of
B
r
was
not restrictive enough, especially near the boundaries of the domain,
and so the rst-order derivative of
B
r
was used instead, i.e.
B
r
.
Formally, these norms, denoted as Q
|
Br |
and Q
|
2
s
Br |
, measure the
weighted mean-square magnitude of
B
r
and the surface Laplacian of
B
r
over the CMB, O, from1960 to 1985, respectively. This leads to
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 45
Q
|
Br |
(x) =
_
1985
1960
_
O
[
B
r
(t. r. . )]
2
dOdt
__
1985
1960
_
O
dOdt
_
1
. (51)
Q
|
2
s
Br |
(x) =
_
1985
1960
_
O
_
2
s
B
r
(t. r. . )
_
2
dOdt
__
1985
1960
_
O
dO dt
_
1
. (52)
where B
r
and its derivatives are computed from the potential of
eq. (5).
4.1.2 Minimal nightside ionospheric E-region currents
It is known that the nightside ionospheric E-region conductivity is
greatly diminished owing to the lack of solar extreme-ultraviolet
(EUV) ionizing radiation, at least at the mid and low latitudes. This
means that the equivalent sheet current density, J
eq
(which is as-
sumed to be equal to the true ionospheric current density) is mini-
mal in these areas. This J
eq
is a toroidal eld residing on the sphere
where the sheet currents ow and may be expressed in terms of a
current or stream function, + (Chapman & Bartels 1940), such that
J
eq
(t. t
mut
. r) = r +(t. t
mut
. r). (53)
where
+(t. t
mut
. r) =
_
1
j
0
2
s=2
4
p=0
p+1
l=p1
|l|+40
k=max (1.|l|)
c
l
ksp
60
n=1
min (n.12)
m=min (n.12)
_
d
lm
kn.e
_
_
2n +1
n +1
_
S
m
nsp.e
(t. t
mut
. r)
_
. (54)
and r is the position vector on the sphere of radius r = a +h.
The theory takes the form of a quadratic norm, Q
Jeq
, which
measures the mean-square magnitude of J
eq
on a spherical sector,
O
s
, xed in dipole MLT longitude (dened as
mlt
=
d
+
p
t
mut
)
over UT. The low-conductivity sector is 8 h wide and centred on
local 1:00. This leads to
Q
Jeq
(x) =
_
T
0
_
Os
J
eq
(
mlt
. t
mut
.
d
)
2
dO
s
dt
mut
_
_
T
0
_
Os
dO
s
dt
mut
_
1
. (55)
The t
mut
integration is facilitated by two assumptions: rst, the time
variation of the F
10.7
index is neglected, rendering J
eq
periodic over
1 yr such that T = 1 yr may be used; and secondly, universal time t
in eq. (9) is treated as t
mut
, which is then used for time integration.
Although the second assumption is best at low and mid latitudes,
it is expected that this will make little difference in the analysis,
especially since Q
|Jeq|
is a soft bound.
The Q
Jeq
norm works in conjunction with the p = 0 terms of
eq. (25) to establish a nightside baseline such that J
eq
is minimized
at those hours. This baseline is a global function, able to adjust to
geographic shifts, which is static on diurnal timescales, but varies
with season. Because there is difculty in separating this functional
behaviour fromthat of main eld SVat satellite altitude, the strength
of the norm is adjusted via the associated
Jeq
such that all p = 0
terms are determined by the norm. It is also for this reason that the
inuence of the norm cannot be greatly reduced in the polar regions
(e.g. via some dipole colatitude weighting function) where J
eq
is
thought to ow at all MLTs.
4.1.3 Smooth ionospheric E-region currents
Recall that in order to resolve the EEJ along the dip equator, QD
degrees of up to k =45 are used. Since the T
l
ksp.e
functions are global,
it is expected that spurious oscillations will be exhibited in the J
eq
morphology. Although the preferred model state for the Q
Jeq
norm
is zero, its inuence is limited in proximity to MLTs of 215 h.
Hence, an additional norm is sought to minimize this roughness on
the dayside, which suggests minimizing the mean-square magnitude
of some function of the second-order horizontal derivatives of J
eq
.
A natural choice is the surface Laplacian of J
eq
. However, this norm
shouldnot interfere withthe p = 0baseline establishedfor J
eq
bythe
Q
Jeq
norm, and so it is restricted to current densities in the p > 0
regime, denoted as J
eq. p>0
. Consequently, the norm may be applied
at all MLTs. Furthermore, it must not interfere with legitimate EEJ
variations near the dip equator nor with ow in the auroral regions.
This may be accomplished by introducing a non-negative weighting
function in
d
, which is smaller in the equatorial and polar regions
and larger at mid-latitudes; for this study sin
8
2
d
is used. Although
a more rigorous approach would use QD colatitude, it is much more
complicated and is left for future work.
Formally, this norm, denoted as Q
2
s
J
eq. p>0
, measures the
weighted mean-square magnitude of the surface Laplacian of J
eq. p>0
on a sphere, O, over UT. This leads to
Q
2
s
J
eq. p>0
(x) =
_
T
0
_
O
_
_
2
s
J
eq. p>0
(t
mut
.
d
.
d
)
_
_
2
sin
8
2
d
dOdt
mut
__
T
0
_
O
sin
8
2
d
dO dt
mut
_
1
. (56)
The same assumptions regarding the t
mut
integration are made here
as well. Note that the surface Laplacian operator multiplies S
m
nsp.e
by
a factor of n(n +1), and so Q
2
s
J
eq. p>0
, which measures
the mean-square magnitude of the deviations from a dipole in MLT
(k >1or l = p) andindependent of D
st
ona sphere at Magsat altitude
(r = a +h
m
with h
m
= 450 km), O, over UT. This leads to
Q
LB
ltd
(x) =
_
T
0
_
O
LB
ltd
(t
mut
.
d
.
d
)
2
dOdt
mut
_
_
T
0
_
O
dO dt
mut
_
1
. (57)
C
2002 RAS, GJI, 151, 3268
46 T. J. Sabaka, N. Olsen and R. A. Langel
where LB
ltd
includes all j
l
ksp
terms in eq. (29) for which k >1 or
l = p.
4.1.5 Minimal ionospheric coupling currents
Recall that the radial component of the meridional coupling currents
of the EEJ are being accounted for in the Magsat observations via
B in eq. (30), requiring QD degrees of up to k = 40. Hence, T
l
ks.j
will be susceptible to instabilities similar to those of T
l
ksp.e
, and
consequently, the associated J
r
for both dawn and dusk will need to
be smoothed. Consider that the inclination of the Magsat orbit was
such that no data were acquired within a cap of half-angle of about
7
2
s
J
eq. p>0
s=0
40
k=1
min (k.4)
l=0
l
ks
60
n=1
min (n.12)
m=min (n.12)
_
d
lm
kn.j
_
n (n +1) S
m
ns.j
(t. r)
_
_
_
. (59)
and r is the position vector on the sphere of radius r =a +h
m
.
Because J
r
has a period of 1 yr, T =1 yr. Again, there are separate
Q
| Jr |
included for both dawn and dusk local times.
4.2 Weighting
In this section the W matrix in eq. (46) will be dened. The com-
ponents of the data noise vector are assumed to reect stationary,
uncorrelated processes, which renders a diagonal W. In particular,
the scalar noise process is considered Gaussian, which is only true
if the process is isotropic.
The data subsets of this study, considered to represent distinct
stochastic populations, are listed as subheaders under the com-
ponent column of Table 3 (see Section 5). Further assumptions
are made concerning the variation of the error processes with re-
spect to orientation: the OHM-1AM and low/mid dipole latitude
(|90
d
| - 50
d
| 50
) OHM-MUL
and Magsat dawn and dusk processes are considered to be isotropic
in the XY-plane; and the Magsat dawn and dusk mid/low dipole
latitude processes are considered to be isotropic in the XZ-plane.
Assignment of the actual variances is accomplished by iteratively
adjusting a starting set of variances similar to s
2
of eq. (48), but cor-
responding to each data subset, until they reach unity. This requires
that a model be tted to the data. Because this is a computationally
intensive venture, a preliminarymodel, denotedA, developedduring
the initial stages of this study is used once for all the data. Adetailed
description of this model is provided by Sabaka et al. (2000). The
a priori data noise variances resulting from this procedure are listed
in the
damping parameter
value. The values selected for
i
are listed in Table 2.
5 RES ULTS
The CM3 model may be examined for self-consistency and for phys-
ical plausibility, the rst of which will now be discussed.
5.1 Residuals and data ts
In this section, the data residuals will be examined. Table 3 lists,
among other things, the unweighted residual statistics, i.e. the mean,
j
r
, and standard deviation about the mean,
r
, of the CM3 model for
each eld measure for each of the data subsets of interest. Also listed
is the a priori standard deviation of the associated data noise,
. j
r
Table 2. Damping parameter values used in the CM3 model.
Norm Damping parameter ()
Q
|
Br |
3.1 10
1
(nT yr
2
)
2
Q
|
2
s
Br |
3.1 10
7
(nT yr
1
km
2
)
2
Q
Jeq
8.4 10
2
(A km
1
)
2
Q
2
s
J
eq. p>0
3.8 10
2
(A km
3
)
2
Q
LB
ltd
6.3 10
4
(nT)
2
Q
| Jr |
dawn 1.3 10
5
(nA m
2
)
2
Q
| Jr |
dusk 1.3 10
5
(nA m
2
)
2
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 47
Table 3. Comparison of unweighted residual statistics (j
r
,
r
, and
in units of nT).
Component CM3 GSFC/CU(12,96) GSFC(8,95-SqM)
Number j
r
r
Number j
r
r
Number j
r
r
OAM
X 4047 0.0 28.5 4048 0.0 26.0
Y 4047 0.0 40.1 4048 0.0 26.0
Z 4047 0.0 33.7 4048 0.0 31.9
OHM-1AM
X 27 094 0.0 15.0 15.0
Y 27 106 0.0 13.1 15.0
Z 27 512 0.0 19.6 15.0
OHM-MUL
|90
d
| - 50
X 54 792 0.0 10.0 11.0 56 963 0.0 11.4 153 544 0.0 14.8
Y 54 846 0.0 9.0 11.0 57 016 0.0 12.1 153 544 0.0 14.2
Z 56 696 0.0 8.5 11.0 55 978 0.0 10.2 153 544 0.0 11.1
|90
d
| 50
X 65 834 0.0 17.0 18.0 65 451 0.0 17.9 74 712 0.0 19.4
Y 65 895 0.0 16.5 18.0 65 487 0.0 17.6 74 712 0.0 17.5
Z 65 693 0.0 23.2 21.0 65 230 0.0 20.8 74 712 0.0 22.4
Magsat corrected dusk
X 11 060 0.3 10.3
Y 11 060 0.3 13.4
Z 11 060 0.3 8.2
B 12 399 0.3 7.5
Magsat dusk
X 9381 0.1 4.6 5.4 9381 0.04 5.3
Y 9321 0.03 5.7 6.8 9321 0.0 6.8
Z 9382 0.1 4.0 5.4 9382 0.1 5.4
B 11 404 0.5 3.6 5.5 11 404 1.04 5.5
Magsat polar dusk
X 7985 1.6 15.8 18.5 7985 2.5 17.6
Y 7988 0.3 16.5 18.5 7988 1.0 17.6
Magsat dawn
X 10 570 0.05 4.3 5.0 10 570 0.3 5.1 10 595 0.4 7.1
Y 10 537 0.0005 4.4 5.4 10 537 0.1 5.2 10 595 0.1 7.4
Z 10 588 2.0 3.4 5.0 10 588 0.8 4.7 10 595 0.2 6.2
B 12 441 0.2 3.6 5.3 12 441 0.6 5.2 12 460 0.6 7.4
Magsat polar dawn
X 8483 0.4 17.4 19.0 8483 1.3 18.3
Y 8445 0.1 18.2 19.0 8445 1.8 19.4
POGO original
B 57 434 2.0 8.0
POGO decimated
B 22 685 0.2 4.2 4.8 22 685 0.0 5.1
POGO pass
B 6754 0.4 5.0 5.8 6754 0.0 6.4
indicate whether a bias exists in the description of the data by the
model. With the exception of Magsat polar dusk Xand Magsat dawn
Z, j
r
have magnitudes well below 1 nT; in fact those of the OHM-
1AMs and OHM-MULs are zero, by virtue of their vector biases.
As expected, the mid/low dipole latitude subsets of Magsat and the
OHM-MULs are tted substantially better than their high-latitude
counterparts. Within the mid/low dipole latitude Magsat dawn and
dusk subsets scalar B is well tted (unaffected by attitude errors) as
is Z (least affected by external contributions) whilst Y is t slightly
worse (possibly caused by attitude error (R. Holme, priv. comm.,
2000) or the result of dynamic variability in the meridional current
system, especially at dusk). The POGO B are tted slightly worse
than Magsat B, which may be caused by their extended distribution
in time and the OHM-1AMs ts are consistent with those of the
bulk OHM-MULs.
The OAMs and OHM-MULs for Tucson are shown in Fig. 5 at two
different timescales, along with the values predicted by CM3. The
top panel shows all available OAM data from 1960 to 1985 along
with the OHM-MULs. The OAMs were predicted using annually
integrated values from the magnetosphere and ionosphere. The pre-
dictions appear to be quite good for all components, especially when
considering that the OAMs were not used when deriving the model.
The bottom panel shows the observed and predicted OHM-MULs
for the quietest day of each month in 1967. The ts are judged to be
good, especially when considering the adjustment for the jump dis-
continuities in X and Z between the quietest days of May and June,
C
2002 RAS, GJI, 151, 3268
48 T. J. Sabaka, N. Olsen and R. A. Langel
Tucson observatory for all years
25000
25100
25200
25300
X
5400
5600
5800
Y
42800
43200
43600
1960 1965 1970 1975 1980 1985
Z
Tucson observatory for 1967
25140
25160
25180
25200
25220
X
5760
5800
5840
Y
43520
43540
43560
43580
43600
Z
J
A
N
F
E
B
M
A
R
A
P
R
M
A
Y
J
U
N
J
U
L
A
U
G
S
E
P
O
C
T
N
O
V
D
E
C
Figure 5. Fits of the CM3 model to the OAM and OHM-MUL vector components of the Tucson observatory. The top panel shows all measured OAMs
(triangles) and OHM-MULs (squares) and their predicted values (grey line), all in nT, over the time span of the model, in years. The bottom panel shows all
OHM-MULs (squares) and the model predictions (black line), all in nT, for the year 1967 (indicated by the box outline in the top panel). The abscissa in the
bottom panel is discontinuous, being comprised of the quietest day of each month over the year, and begins at 00 UT for each day.
which are caused by a change in external eld strength (ring-current
level).
Similar plots are provided in Fig. 6 for the Huancayo observa-
tory, which is located under the EEJ. The top panel shows a similar
behaviour as in the case of Tucson, except that the Z predictions
are excessively smooth, especially near the end points. This may be
caused by the smoothing of the main eld SV in Z, which is dom-
inating the small-Z component at this dip equatorial station. The
bottom panel is now for the year 1966 and shows good ts to all
three components. Note the adjustment in X for the jump disconti-
nuity between the quietest days of June and July caused by different
levels of ring-current activity.
While the Tucson and Huancayo plots show that the model sat-
isfactorily predicts the analysed OHM-MUL data in conned en-
velopes and the integrated effects of OAMs across the domain,
it is also important that the model be able to extrapolate its full
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 49
Huancayo observatory for all years
27000
27300
27600
27900
28200
X
1200
1500
1800
2100
2400 Y
900
950
1000
1050
1100
1960 1965 1970 1975 1980 1985
Z
Huancayo observatory for 1966
27920
28000
28080
X
2160
2200
2240
Y
1040
1060
1080
Z
J
A
N
F
E
B
M
A
R
A
P
R
M
A
Y
J
U
N
J
U
L
A
U
G
S
E
P
O
C
T
N
O
V
D
E
C
Figure 6. Fits of the CM3 model to the OAM and OHM-MUL vector components of the Huancayo observatory. The top panel shows all measured OAMs
(triangles) and OHM-MULs (squares) and their predicted values (grey line), all in nT, over the time span of the model, in years. The bottom panel shows all
OHM-MULs (squares) and the model predictions (black line), all in nT, for the year 1966 (indicated by the box outline in the top panel). The abscissa in the
bottom panel is discontinuous, being comprised of the quietest day of each month over the year, and begins at 00 UT for each day.
resolving power to poor control areas. To this end, the rms ts to
the OHM-MULs from three observatories, Addis Ababa (equato-
rial), Gnangara (mid-latitude, southern hemisphere) and Thule (high
latitude), were computed in 1 yr bins from 1960 to 1985. The re-
sults are shown in Fig. 7, with the rms X, Y and Z values indicated
by squares, triangles and circles, respectively. The data segments
within the POGO and Magsat coverage envelopes were analysed
in the model and these are indicated with a shaded backdrop. All
values ramp up near 1960 and most ramp up slightly near 1985.
However, in the satellite gap between 1972 and 1979, only the Y
component of Addis Ababs shows any excursion more severe than
at any other time. In fact, it could be argued that several compo-
nents are actually better tted in this gap than during the satellite
coverage envelopes. Note that Addis Ababa X and Thule Z show
the highest variance, as would be expected from their geographic
positions.
C
2002 RAS, GJI, 151, 3268
50 T. J. Sabaka, N. Olsen and R. A. Langel
Addis Ababa
0
20
40
60
80
1960 1965 1970 1975 1980 1985
0
20
40
60
80
1960 1965 1970 1975 1980 1985
0
20
40
60
80
1960 1965 1970 1975 1980 1985
Gnangara
0
20
40
60
80
1960 1965 1970 1975 1980 1985
0
20
40
60
80
1960 1965 1970 1975 1980 1985
0
20
40
60
80
1960 1965 1970 1975 1980 1985
Thule
0
20
40
60
80
1960 1965 1970 1975 1980 1985
0
20
40
60
80
1960 1965 1970 1975 1980 1985
0
20
40
60
80
1960 1965 1970 1975 1980 1985
Figure 7. The rms ts to OHM-MULs from the Addis Ababa (equatorial), Gnangara (mid-latitude, southern hemisphere) and Thule (high-latitude) obser-
vatories, computed in 1 yr bins from 1960 to 1985. The rms X, Y and Z values are indicated by squares, triangles and circles, respectively. The shaded areas
delineate the POGO and Magsat coverage envelopes where the OHM-MULs were actually analysed in the model. All ordinates are in nT.
The top panel of Fig. 8 shows the t to B for a particular pass
of POGO data that was included in the model. This pass crosses
the geographic equator at 59
geographic
latitude. The coupling currents member of the Ysuite shows signif-
icant toroidal eld presence at lowlatitude, which is caused by radial
meridional currents connected with the EEJ impinging the Magsat
sampling shell. This model (probably the rst to include the global,
non-potential contribution at satellite altitude) is able to t this
feature.
5.2 Parameter uncertainty and separability
As with any inverse problem, one must address the issues of pa-
rameter uncertainty and separability. These ultimately reect upon
the quality and distribution of data and any additional regulariz-
ing information from which the parameters are inferred. It must
be stressed that in this study the measures of uncertainty and sep-
arability are statistical quantities derived from a linear covariance
analysis, which is applicable here since the non-linearities owing to
scalar data are very slight and noise distributions are assumed to be
Gaussian. The a posteriori parameter error-variances may be taken
directly from the diagonals of the calibrated error-covariance ma-
trix, C
x
, introduced in Section 4. However, if signicant off-diagonal
elements exist, then the true resolving power of the data may be ob-
fuscated and will most likely lie along linear combinations of the
parameters. A principal component analysis, or eigenvalue decom-
position, of C
x
will, in fact, reveal these combinations and give a
much clearer viewof what is known or unknown by the data. Fig. 10
shows the square-root spectrum (principal standard deviations) of
C
x
, which may be read directly as the uncertainties (in nT) of the
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 51
POGO B residual suite
-60
-40
-20
0
20
40
Magnetosphere
-60
-40
-20
0
20
40
Ionosphere
-60
-40
-20
0
20
40
-90 -60 -30 0 30 60 90
Lithosphere
Magsat Z residual suite
-40
-20
0
20
40
Magnetosphere
-40
-20
0
20
40
Ionosphere
-40
-20
0
20
40
-90 -60 -30 0 30 60 90
Lithosphere
Figure 8. Fit of the CM3 model to the scalar (B) values of a particular POGO pass used to derive the model and to the Z component of the Magsat dusk pass
263, which was not used in deriving the model. The POGO pass crosses the geographic equator at 59
W. The
top panel is a suite of residual plots for the POGO B data as a function of latitude. The progression is from the top to bottom member, with a given member
showing residuals, LB (squares), with respect to the main eld plus all preceding labelled elds, B
cum
, and the component (black line) of the predicted, currently
labelled eld, B
new
, in the direction of B
cum
. The bottom panel is a similar suite of residual plots for the Magsat Z data, except the black line now represents
the Z component of B
new
. All ordinates are in nT.
corresponding eigenvectors, which represent these linear combina-
tions. Arrows indicate boundaries in which OHM bias (Biases),
coupling-current (CCs), core/lithospheric (C/L), magnetospheric
(Mag), or ionospheric (Ion) parameter directions account for at least
1 per cent of some of the eigenvector lengths, and so, give a sense as
to how well these groups are known. First, the shaded area indicates
the region in which 99.99 per cent of the total model error-variance is
concentrated; it encompasses 1189 DOFs, and combinations outside
of this region are known to better than 0.77 nT. This region is dom-
inated almost exclusively by uncertainties in OHM biases. In par-
ticular, the least-known combinations of parameters involve several
Antarctic stations, with the South Pole Z biases having the highest
uncertainty at 225 nT. There is a subregion in which eigenvectors
span both biases and coupling-current parameters, indicating the
presence of a slight correlation between these two non-potential
elds. The eigenvectors spanning the parameters of the remaining
source elds are all known very well. However, this does not mean
that the individual parameters themselves are well separated. In fact,
this would translate to perfect separability only if the eigenvalues
(principal variances) were equal for all of these spanning eigenvec-
tors; Fig. 10 clearly shows that this is not the case. One must look
to the correlation matrix associated with C
x
to address separability.
C
2002 RAS, GJI, 151, 3268
52 T. J. Sabaka, N. Olsen and R. A. Langel
Magsat X residual suite
-40
-20
0
20
Magnetosphere
-40
-20
0
20
Ionosphere
-40
-20
0
20
Coupling currents
-40
-20
0
20
-90 -60 -30 0 30 60 90
Lithosphere
Magsat Y residual suite
-20
0
20
40
Magnetosphere
-20
0
20
40
Ionosphere
-20
0
20
40
Coupling currents
-20
0
20
40
-90 -60 -30 0 30 60 90
Lithosphere
Figure 9. Fits of the CM3 model to the X and Y components of the Magsat dusk pass 263. The top and bottom panels are suites of residual plots for the X
and Y data, respectively, from this pass as a function of latitude. The progression is from the top to bottom member, with a given member showing residuals
(squares) with respect to the main eld plus all preceding labelled elds and the X and Y components (black lines) of the predicted currently labelled eld. Note
that coupling currents are now included. All ordinates are in nT.
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 53
Figure 10. Principal standard deviations fromthe calibrated error-covariance matrix, C
x
, of CM3. The shaded area indicates the region in which 99.99 per cent
of the total model error-variance is concentrated. Arrows indicate boundaries in which OHM bias (Biases), coupling-current (CCs), core/lithospheric (C/L),
magnetospheric (Mag) or ionospheric (Ion) parameter directions account for at least 1 per cent of some of the eigenvector lengths.
If the (i j )th element of the parameter error-covariance matrix is
given by
(C
x
)
i j
=
a
( x
i
x
i
)
a
( x
j
x
j
). (60)
where
a
( x
i
x
i
) is the gradient of the ith parameter error, x
i
x
i
,
with respect to the data model, a, then the (i j )th element of the
parameter correlation matrix, R
x
, is given by
(R
x
)
i j
=
a
( x
i
x
i
)
a
( x
j
x
j
)
a
( x
i
x
i
)
a
( x
j
x
j
)
= cos
i j
. (61)
where
i j
is the angle between the gradients of x
i
x
i
and x
j
x
j
. In this study, parameter separability is dened as the degree of
independence of one parameter error from another with respect to
data perturbations. Indeed, if two parameters are perfectly correlated
or anticorrelated, then there does not exist a data perturbation that
affects one but not the other. This leads to a singular C
x
matrix.
Conversely, if two parameters enjoy some degree of separation, then
there exists a data perturbation, which affects one to some degree but
not the other. Thus, the correlation matrix gives a direct statistical
measure of this degree of separability.
Elements of the CM3 correlation matrix for which the absolute
values exceeded 0.7 were examined for evidence of parameter sep-
arability problems. This threshold was chosen since it corresponds
to
i j
45
d
| - 50
d
| 50
Br |
25 839.3 2535 2140.6 394.4 65.51
Q
|
2
s
Br |
6795.1 2535 45.3 2489.7 2.73
Q
Jeq
11 108.1 5520 1591.4 3928.6 2.83
Q
2
s
J
eq. p>0
Br |
, Q
|
2
s
Br |
, Q
Jeq
, Q
2
s
J
eq. p>0
,
and Q
LB
ltd
, respectively, while Q
| Jr |
dawn and dusk both resolve
about 3 per cent of the parameters.
The data subset and prior mist factors are listed in the last col-
umn. These are ratios of the observed to the expected errors of
the model, and hence, give an indication of the relative inuence
of the information associated with each factor. These factors are
not independent and the adjustment of one will affect the others
in complicated ways (see Seber & Wild 1989) for discussions on
the related topic of the iteratively reweighted least-squares (IRLS)
method). However, values above (below) unity would suggest a rela-
tive decrease (increase) in inuence of that information is warranted
to achieve calibrated error bounds. The most noteworthy deviations
from unity come from the Q
|
Br |
and Q
LB
ltd
accounts for less than 1 per cent of the total variance, its inuence
is considered minor. As for Q
|
Br |
, a reduction in its inuence would
indeed allowfor much better ts to the OHMs used in the model, but
the resulting ts to the OAMs would be unrealistically poor. Finally,
the total data mist (s
2
) is 0.87. This suggests that the overall data
uncertainties be reduced by about 7 per cent to achieve calibration.
Indeed, this is an efcient eld model when considering its data
mist in light of its data-to-parameter ratio of 36 compared with
that of the GUFM1 model of Jackson et al. (2000) (a mist of 1.16
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 57
and a data-to-parameter ratio of 10) and the OIFM model of Olsen
et al. (2000) (a mist of 1.01 and a data-to-parameter ratio of 34).
6 DI S CUS S I ON
In this section the examination of the CM3 model moves to more
physical grounds. The salient features of the model are compared
with, and assessed in light of, other works and with the known
physics of the near-Earth magnetic eld.
6.1 Previous comprehensive models
Before a discussion of the various source elds is undertaken, it is
instructive to compare the unweighted residual statistics of the CM3
model with its predecessor GSFC(8,95-SqM), introduced in Sec-
tion 1.3 and for which details are given in Sabaka et al. (2000) and
in the original paper. The statistics for this model, as well as for the
GSFC/CU(12,96) model introduced in Section 4.3, are given in Ta-
ble 3. Note that only the satellite data sets are identical between CM3
and GSFC/CU(12,96), otherwise the statistics are arranged by qual-
itatively similar groupings. With the exception of the high-latitude
Z component for CM3, the
r
for all components of the OHMs
generally decrease with increasing model sophistication. The CM3
Z degradation is caused by SV smoothing in B
r
. This is slightly
misleading, however, since signicant OHM-MUL data discrep-
ancies were detected and corrected between the GSFC/CU(12,96)
and CM3 analyses (note the differences in data counts). Since errors
in the OAM data were initially suspected, these were replaced by
OHM-1AMvalues. This also expedited the analysis since OAMs are
difcult to analyse properly if the integrated nature of their measure-
ments are to be treated properly. The results are that the CM3 ts for
OHM-1AM and OHM-MUL are very consistent as well as far su-
perior to the OAMts of GSFC(8,95-SqM) and GSFC/CU(12,96).
The CM3 ts to the OAMs will be discussed shortly. The low-/mid
-latitude Magsat data sets also exhibit a decrease in
r
as the mod-
els progress, even when considering the corrected versus uncor-
rected dusk data sets. For CM3,
r
for the vector components of the
mid-/low-latitude Magsat data are actually less than the 6 nT ux-
gate accuracy quoted by Langel & Hinze (1998). However, a more
modern treatment, which takes into account such things as attitude
anisotropy (see Holme & Bloxham 1996), would show that this is
probably an overestimate of the actual error. Therefore, tting better
than 6 nT does not imply an overtting of the data. The
r
for the
POGO data sets also decrease dramatically from GSFC/CU(12,96)
to CM3. This is because the SV inuences POGO much more than
Magsat and the SV is much improved because of the OHM data er-
ror corrections. As expected, the CM3 model is providing superior
data ts to those of its predecessors.
6.2 Main eld secular variation
It is essential that the CM3 SV model be plausible not only because
the disjoint Magsat and POGO data envelopes must be correctly
linked for proper parameter determination, but because the entire
time domain of the model is of scientic interest. To assess the SV,
the model will be compared with the GUFM1 model of Jackson
et al. (2000) in terms of Gauss coefcients (eq. 5) and ts to OAMs,
which were not used when deriving the CM3 model. The GUFM1
model, the successor to the UFM1 and UFM2 models of Bloxham
& Jackson (1992), spans four centuries (15901990) and is derived
from a vast historical data compilation. Because the main focus of
this model is describing the SV at the CMB, it is well suited as a
standard for comparison. The model considers internal eld sources
only and represents the Gauss coefcients as B-spline expansions
in time. In Fig. 14, the rst six Gauss coefcients from CM3 (solid)
and GUFM1 (dashed) are shown over the time span of CM3. The
agreement is quite good, but deviations do occur mostly near the
boundaries of the interval. Fig. 15 shows the rates of change of
these same coefcients. Both models show the same SV to rst
order, however, there are some signicant differences. Recall that
GUFM1 includes no external eld contributions nor their induced
effects. It is therefore of interest to see whether at least some of
the differences may be attributable to aliased long-period induction
effects in GUFM1. Of course, this is as much a test for CM3 in its
ability to separate core and external time variations. The axial dipole
term, in dipole coordinates, from CM3 was subtracted from that of
GUFM1 through time and is shown in Fig. 16. It was felt that this
component might reveal an aliased signal if one existed. While one
could argue from visual inspection that power exists at about the 20
and perhaps even the 10 yr period, a comparison with the predicted
signal from the induced portion of CM3 was unconvincing in ac-
counting for the bulk of this difference. Furthermore, this difference
is no doubt sensitive to the amount of damping applied to the SV.
The quality of the SV model may also be addressed by its ability
to predict long time-series of observatory data, such as the OAMs,
which were not analysed in the model. Table 5 compares the rms
predictions with 4047 OAM vector observations from GUFM1 and
CM3. An additional column shows the effects of not including exter-
nal eld contributions from CM3. The rms values are indeed very
close, but the full CM3 model does better for the X and Z com-
ponents. Note that external elds, particularly the magnetospheric
eld, provide the greatest improvement in the Xcomponent for CM3.
Additional measures of the CM3 SVmodel are provided in Table 6.
6.3 Core at 1980 and lithospheric elds
The R
n
spectrum of Lowes (1974) and Mauersberger (1956), which
measures the mean-squared magnitude of the internal eld over a
sphere at a particular epoch per spherical harmonic degree, provides
useful information concerning global properties of the spatial char-
acter of the model. The spectrumfor CM3, computed at r = 6371.2
km and epoch 1980, shows the expected distinct change in slope
of ln (R
n
) around n =14, where the core gives way to crustal dom-
ination, as noted by Langel & Estes (1982). However, there also
appears to be a distinct noise oor emerging prior to n = 50. This
type of behaviour is also reported in the spectrum of the M07AV6
model of Cain et al. (1989b). This suggests at least three spec-
tral regimes: low-degree core-dominated (S
1
), mid-degree crustal-
dominated (S
2
) and high-degree noise-dominated (S
3
). Since the as-
signment of boundaries between the regimes can be somewhat sub-
jective, especially when done visually, a more objective approach
was developed in which a three-segment best-tting linear piece-
wise regression (BFLPR) to unweighted ln (R
n
) (excluding ln (R
1
))
was performed as a function of degree partitioning, i.e. the degree
boundaries were chosen that minimized the total mist. Note that a
segment fromone regime was not extrapolatedintothe other regimes
during the tting procedure; a process that does not reect the true
spectral overlap of the sources. However, for purposes of compari-
son this will sufce. The resulting BFLPR is
ln (R
n
) =
_
_
(1.27 0.07)n +(20.8 0.7) for S
1
.
(0.05 0.02)n +(1.8 0.7) for S
2
.
(0.11 0.03)n (1.0 1.7) for S
3
.
(66)
C
2002 RAS, GJI, 151, 3268
58 T. J. Sabaka, N. Olsen and R. A. Langel
Main field: CM3 (solid), GUFM1 (dashed)
-30400
-30350
-30300
-30250
-30200
-30150
-30100
-30050
-30000
-29950
-29900
g
1
0
1960 1965 1970 1975 1980 1985
-2150
-2100
-2050
-2000
-1950
g
1
1
1960 1965 1970 1975 1980 1985
5500
5550
5600
5650
5700
5750
h
1
1
1960 1965 1970 1975 1980 1985
-2050
-2000
-1950
-1900
-1850
-1800
-1750
-1700
-1650
-1600
-1550
g
2
0
1960 1965 1970 1975 1980 1985
2990
2995
3000
3005
3010
3015
3020
3025
3030
3035
3040
g
2
1
1960 1965 1970 1975 1980 1985
-2200
-2150
-2100
-2050
-2000
h
2
1
1960 1965 1970 1975 1980 1985
Figure 14. First six Gauss coefcients from CM3 (solid) and the GUFM1 (dashed) model of Jackson et al. (2000), all in nT, over the time span of CM3.
where S
1
, S
2
and S
3
are found to correspond to degree ranges 214,
1542 and 4365, respectively. This t is shown, along with the raw
R
n
values, in the top panel of Fig. 17. Points of intersection between
regressions for S
1
and S
2
, S
1
and S
3
, and S
2
and S
3
occur at n = 14.4,
15.7 and 43.1, respectively, as compared with the S
1
and S
2
value
of n = 14.2 given by Cain et al. (1989b). Under the assumptions
of optimal (Wiener) ltering (Press et al. 1992), the noise spectrum
is considered uncorrelated with that of the core and the crust and is
hypothesizedtofollowthe same trendat lower degrees as established
at higher degrees. The bottom panel of Fig. 17 shows the BFLPR to
ln (R
n
), over the same core- and crustal-dominated regimes of the
top panel, after subtraction of R
n
noise values extrapolated from the
linear regression for the high-degree spectrum. These regressions
are given by
ln (R
n
) =
_
(1.27 0.07)n +(20.8 0.7) for S
1
.
(0.010 0.023)n +(2.7 0.7) for S
2
.
(67)
While the noise correction has had a negligible affect on the core-
dominated spectrum, it has removed practically all of the slope in
the linear regression for the crustal-dominated portion, rendering it
nearly level. The radii, R
= R
d. (68)
where ln (d) is the slope of the linear regression. Applying this to the
corrected core spectrum yields a levelling depth of 116 125 km
below the seismic coremantle boundary at r =3485 km. This falls
generally in between the values of 174 km given by Langel & Estes
(1982) and 8046 kmgiven by Cain et al. (1989b), although the un-
weighted error envelope encompasses these values as well as several
kilometres above the CMB. The levelling depth for the crustal spec-
trum is 31 74 km below the Earths surface, where the Cain et al.
(1989b) value is 21 km. This exercise, however, should not be taken
out of context given the overly simplifying assumptions that have
been made. There is evidence, for instance, fromthe statistical mod-
els of crustal magnetization of Jackson (1994) that the crustal power
may be increasing with degree over this range and only at higher
degrees (n = 200) does it begin to fall-off. In addition, as mentioned
earlier, spectral overlap was not considered when establishing the
regimes. However, when the overlap is properly accounted for, the
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 59
Secular Variation: CM3 (solid), GUFM1 (dashed)
18
20
22
24
t
g
1
0
1960 1965 1970 1975 1980 1985
4
6
8
10
12
t
g
1
1
1960 1965 1970 1975 1980 1985
-24
-22
-20
-18
-16
-14
-12
-10
-8
-6
-4
-2
0
t
h
1
1
1960 1965 1970 1975 1980 1985
-26
-24
-22
-20
-18
-16
-14
t
g
2
0
1960 1965 1970 1975 1980 1985
0
2
4
t
g
2
1
1960 1965 1970 1975 1980 1985
-16
-14
-12
-10
-8
-6
-4
-2
t
h
2
1
1960 1965 1970 1975 1980 1985
Figure 15. Rate of change of rst six Gauss coefcients from CM3 (solid) and the GUFM1 (dashed) model of Jackson et al. (2000), all in nT yr
1
, over the
time span of CM3.
-6
-4
-2
0
2
4
6
8
g
(
n
T
)
1960 1965 1970 1975 1980 1985
time (yrs)
Figure 16. Axial dipole term (g
0
1
), in dipole coordinates, of CM3 subtracted from that of GUFM1.
levelling depth for the crustal spectrumbecomes 122 37 kmabove
the Earths surface; a value close to the ionospheric E-region sheet
current height. For comparative purposes, however, these exercises
show that the R
n
spectrum of CM3 is in good agreement with previ-
ous work, but see Walker & Backus (1997) for a more sophisticated
treatment.
In light of the discussion on R
n
, the lithospherically dominated
portionof the internal eldmodel is takentobe n =1542andthe as-
sociated B
r
map at satellite altitude is shown in Fig. 18. The LZ map
for n =1565 of Ravat et al. (1995) (referred to as the RLPAAmap)
will be usedfor comparison. It is derivedfromMagsat data usingvar-
ious data processing techniques as well as techniques for modelling
C
2002 RAS, GJI, 151, 3268
60 T. J. Sabaka, N. Olsen and R. A. Langel
Table 5. Comparison of rms predictions to OAM data (units in nT).
Component Number rms
GUFM1 CM3
All No external
X 4047 17.71 17.48 18.09
Y 4047 21.27 21.45 21.47
Z 4047 24.55 24.49 24.53
Table 6. CM3 secular variation measures at the CMB from 1960 to 1985.
rms radial secular variation 1220 nT yr
1
rms radial secular acceleration 28 nT yr
2
rms surface Laplacian of radial 1.2 10
2
nT yr
1
km
2
secular variation
and removing ionospheric elds, and a covariant spherical harmonic
analysis procedure to isolate common dawn and dusk lithospheric
anomaly features. First, both maps are in excellent agreement with
regard to the shapes and locations of the major well-known anoma-
lies such as Kiruna (A), Kursk (B), Bangui (C), Gulf of Mexico (D),
Kentucky (E), Australian Bight (F), etc. (letters identify the anomaly
locations in Fig. 18). In general, however, the anomaly magnitudes
R
n
spectra for CM3 at Earths surface
10
0
10
1
10
1
10
2
10
3
10
4
10
5
10
6
10
7
10
8
10
9
10
10
n
T
2
0 10 20 30 40 50 60 70
Degree
Uncorrected
10
0
10
1
10
1
10
2
10
3
10
4
10
5
10
6
10
7
10
8
10
9
10
10
n
T
2
0 10 20 30 40 50 60 70
Degree
Corrected
Figure 17. A comparison of R
n
spectra for the CM3 model at r = 6371.2 km for epoch 1980 corrected (bottom) and uncorrected (top) for high-degree
noise contamination. The solid lines in the uncorrected plot are the three best-tting linear piecewise regressions (BFLPR) to ln (R
n
) as a function of degree
partitioning. The solid lines in the corrected plot are the BFLPR, over the same low- and mid-degree segments, to ln (R
n
) after subtraction of values extrapolated
from the linear regression for the high-degree segment in the uncorrected spectrum.
from CM3 appear larger, as much as 3050 per cent in some cases.
This may be an affect of applying Kaiser ltering (Kaiser 1974) with
a cut-off of 12 000 km to the equatorial/mid-latitude Magsat vector
data used in deriving RLPAA, resulting in diminished amplitudes.
This type of along-track ltering may also remove northsouth in-
formation from Magsat passes. Some evidence of this may be seen,
for instance, when comparing the maps in the eastern Australian
basin region (G) (see Purucker et al. 1997) and the region of the
IzuBonin subduction zone (H) (30
N, 140
N and below 60
N or below about 83
S, 109
S (bottom) of the
CM3 lithospheric contribution (degrees 1542) to B
r
on the sphere r =
6771.2 km (stereographic projections). Letters identify particular anomalies
referenced in the discussion.
polar maps also reveal quite a lot of structure in the CM3 litho-
spheric model that lies near the geographic poles, outside the region
of Magsat data coverage. Recall that CM3 includes data from the
POGO satellites OGO-2, OGO-4 and OGO-6, for which the orbit
inclinations were 87.3
, 86.0
and 82.0
half-
angles. Abetter standard for comparison in these regions is therefore
the POGO-derived anomaly maps of Langel (1990). Since these are
maps of scalar anomalies reduced to pole, they should be very close
to LZ and LB
r
near the north and south poles, respectively. The
strong positive-negative anomaly pair entwined at the north pole (K)
also appears on the POGO map, except here the intensities are
stronger than can be accounted for by upward continuation to the
500 km level of the POGO maps. Near the south pole there is a
strong negative lineation (L) that parallels the 135
315
merid-
ian and corresponds very well with the Trans-Antarctic mountain
chain. This feature is diffuse and fragmented at best in the POGO
map. Again, the strong negative anomaly over the south magnetic
dipole (J) is missing in the POGOmap. Finally, there is a strong new
feature (M) in the CM3 map located at about 86
S, 90
E with a peak
value of almost 19 nT. While it has no counterpart in the POGOmap,
the CM3 model does include POGO measurements over a portion
of this region. However, since the POGOdata gap also transects this
feature, an interpretation must be suspended until further analysis
is performed.
6.4 Ionospheric eld
The most dominant of the expansion coefcients describing the
ionospheric E-region currents are { c
0
1.0.0
}, { c
1
2.0.1
} and { c
2
3.0.2
}.
At sunspot minimum, where typical values of F
10.7
= 66.0 10
22
W m
2
Hz
1
, these have values of { c
0
1.0.0
} = 4.0 nT, { c
1
2.0.1
} =
C
2002 RAS, GJI, 151, 3268
62 T. J. Sabaka, N. Olsen and R. A. Langel
5.6 nT and { c
2
3.0.2
} = 3.2 nT, while at sunspot maximum the
coefcient magnitudes are about twice as large. All are local time
terms (l = p) with no seasonal variation (s = 0), the maximum
amplitudes of which are centred at the subsolar point (real parts) and
contribute to the two major Sq foci of antipolarity in the northern and
southern hemispheres (k l = 1). The rst coefcient describes the
general polarity between the north and south foci at all local times
( p =0), the second describes the regular diurnal variation ( p =1)
of such a eld through the course of a solar day, while the third
imparts some semi-diurnal variation ( p = 2) to this.
The dominance of { c
1
2.0.1
} can be seen in global maps of the
E-region equivalent current function + of eq. (54). In Fig. 20, +
is shown at MLT noon on 1980 March 21, but for different values
of MUT. A value of F
10.7
=140.0 10
22
W m
2
Hz
1
, an average
over the time span of this model, was used to generate the maps.
The dual Sq foci are indeed the major features, showing a slightly
asymmetric current load (rather than the symmetric one expected
for equinox) owing in oppositely directed vortices in the north-
ern (counterclockwise) and southern (clockwise) hemispheres in
accordance with eq. (53). The total current owing in the northern
(southern) vortex is the same (slightly lower) than that reported by
Malin & Gupta (1977), who also included p = 0 terms. However,
their data were acquired during IGYwhen F
10.7
levels were elevated
and their analysis was unable to distinguish between ionospheric and
magnetospheric effects. Besides these, there are several other items
of interest: a signicant decrease in + exists for much of the
darkside hemisphere at all MUT, indicating that the Q
Jeq
constraint
is effective; the boundary between the two foci is coincident with
the dip equator at all MUT, i.e. current ows tangential to the dip
equator at and near local noon, thus afrming the utility of the QD
constraints; there is also some amplitude and shape modulation with
00 MUT
TUC
HUA
06 MUT
TUC
HUA
12 MUT
TUC
HUA
18 MUT
TUC
HUA
Figure 20. Global maps of ionospheric E-region equivalent current function + for 1980 March 21 (Mollweide projection). A value of F
10.7
=140.0
10
22
Wm
2
Hz
1
, an average over the time span of this model, was used to generate the maps. Recall from eq. (54) that + is dened on the sphere r =
6481.2 km. Each of the four panels is centred on noon MLT, but for different MUTs. The associated induced contribution is not included here. A 20 kA current
ows between the contours. Locations of the Tucson (TUC) and Huancayo (HUA) observatories are shown.
MUTthat is beyond what is inherent in the QDconstraints and which
is attributable to non-local time variation; and there is a marked in-
crease in + parallel to the dip equator at and near local noon
for all MUT, revealing an enhanced eastward current ow, which is
in fact the EEJ.
In Fig. 21, + is shown at noon MLT in the centre of the gure,
but for the equinoxes and solstices. A distinct seasonal variation
is evident, as expected. Very clearly, the summer focus occurs at
an earlier local time than the winter focusa known feature that
persists at all MUTs. Though not apparent in Fig. 21 at noon MUT,
similar plots at midnight MUTshowa minimumnorthern (southern)
focal current intensity in December (June) and a maximumnorthern
focal intensity in June, as expected, but a maximum southern focal
intensity in September. This general peculiarity in focal intensity
variation with season is probably attributable to the Q
Jeq
constraint
that imposes circuit closure on the dayside hemisphere and, together
with the static ionospheric eld, may very well be modifying the
expected variation. Seasonal oscillations in + focal intensities with
respect to a preliminary CM3 model including no Q
Jeq
constraint
are indeed of higher amplitude, conrming this idea.
These E-region current signatures are also clearly seen in the eld
intensity of the POGO pass in Fig. 8 and the eld components of
Magsat dusk pass 263 in Figs 8 and 9. Note the accentuation of the B
minimumat noon over the dip equator as well as for X fromthe dusk
extent of the EEJ. The OHM component plots from Tucson (Fig. 5)
and Huancayo (Fig. 6) also reveal something of the daily, seasonal
and spatial structure of the E-region currents. Noting that the daily
segments begin at 00 UT, it can be seen that Tucson X and Z both
have a single negative spike at local noon, while Y spikes positive
just before and negative just after local noon. This agrees with the
location of Tucson being below and just north of the northern focus
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 63
MAR 21
TUC
HUA
JUN 21
TUC
HUA
SEP 21
TUC
HUA
DEC 21
TUC
HUA
Figure 21. Global maps of ionospheric E-region equivalent current function + dened on the sphere r =6481.2 km (Mollweide projection). A value of
F
10.7
= 140.0 10
22
W m
2
Hz
1
, an average over the time span of this model, was used to generate the maps. Each of the four panels is centred on noon
MLT and MUT, but for different seasons, i.e. 1979 December 21 and March 21, June 21 and 1980 September 21. The associated induced contribution is not
included here. A 20 kA current ows between the contours. Locations of the Tucson (TUC) and Huancayo (HUA) observatories are shown.
of + in Fig. 20. At Huancayo, X spikes strongly positive (the EEJ)
and Z spikes negative at local noon, which agrees with its location
below and perhaps slightly north of the dip equator as shown in
Fig. 20. In addition, Tucson Y and Z and all Huancayo components
indicate an intensication of the ionospheric contribution to the
diurnal signal during their respective summers.
6.5 Magnetospheric eld
As expected, the static termalong the dipole axis, {j
0
1.0.0
}, is by far
the most dominant, having a value of 21.4 nT. The magnitudes of the
annual and semi-annual variations along this axis are about 13 and 7
per cent of this magnitude, respectively. Although the mean tilt of the
non-D
st
-dependent magnetospheric dipole eld over season is only
about 2.7
at the
solstices, i.e. |{j
1
1.1.1
} +{j
1
1.1.1
}| = 4.8 nT.
As for higher multipole terms, some have been detected in previ-
ous studies. Malin & Mete Isikara (1976) studied the annual varia-
tion of midnight values of observatory data and found contributions
to j
0
1.1.0.g
, j
0
2.1.0.g
and j
1
2.1.1.g
(where the subscripted g indicates
coefcients in geographic spherical-polar coordinates) that they at-
tributed to seasonal movement of the ring-current relative to the
equatorial plane. Note that their analysis is unable to distinguish
annual variation of the magnetospheric part (j
m
nsp.g
) from that of the
ionospheric part (c
m
nsp.g
). If the real part of a particular coefcient,
say {c
m
np
(t )}, is Fourier analysed for annual periodicity such that
_
c
m
np
(t )
_
= A cos(
s
t ). (69)
then Table 7 summarizes their results under the M&MI column.
Comparable results are also shown for CM3 and a preliminary CM3-
type model without Q
LB
ltd
) with very
tight agreement for {j
0
2.0.g
+c
0
2.0.g
} and less so for {j
0
1.0.g
+c
0
1.0.g
}.
The amplitudes, A, of CNXS are about half those of M&MI, while
those of CM3 are even less. This decreased amplitude in the annual
variation of the CM3-type models is probably due in part to the
Q
Jeq
smoothing that operates on ionospheric seasonal variations
at all p values. However, since the magnetospheric annual variations
are not compensatingfor this smoothing, thenperhaps the bulkof the
M&MI amplitudes is actually of ionospheric origin. The apparent
degradation fromCNXS to CM3 with respect to the M&MI ndings
is probably the result of the Q
LB
ltd
1.6 231
1.9 220
{j
0
2.0.g
+c
0
2.0.g
} 6.0 352
2.8 349
3.6 350
{j
1
2.1.g
+c
1
2.1.g
} 1.3 20
0.2 18
0.4 31
{j
0
1.0.g
} 2.7 222
3.1 220
{j
0
2.0.g
} 0.1 353
1.4 343
{j
1
2.1.g
} 0.04 81
0.4 59
C
2002 RAS, GJI, 151, 3268
64 T. J. Sabaka, N. Olsen and R. A. Langel
Table 8. Magnetospheric daily variation (nT) from T87We Peredo et al. (1993) compared with CM3-type models.
T87We (K
p
= 0) T87We (K
p
= 3) CM3 CNXS
{j
m
n
(t
s
. t
mut
)} Season A A A A
Diurnal {j
1
1
} December 3.0 180
3.3 180
4.0 185
3.7 182
Equinox 0 0 0.7 55
2.1 162
{j
1
2
} December 1.2 359
2.3 359
0.2 343
1.3 356
2.7 357
0.3 348
3.2 0
{j
0
1
} December 1.3 181
1.9 181
0.1 329
0.6 218
Equinox 0 0 0.2 97
0.5 58
Semi-diurnal {j
2
2
} December 0.04 359
0.1 359
0.1 346
0.3 328
2.6 357
{j
1
1
} December 0.8 181
1.1 181
0.1 171
0.7 177
0.7 174
0.1 230
0.6 151
_
j
m
n
(t
s
. t
mut
)
_
= A cos (
p
pt
mut
). (70)
then Table 8 gives a comparison of the amplitude, A, and phase,
, of some coefcients of interest from T87We, CM3 and CNXS.
The agreement of the amplitudes and phases of {j
1
1
} and {j
1
2
} in
December and at equinox between the T87We and CNXS models
is considered good, except in the case of the diurnal variation of
the former coefcient at equinox, which is negligible in the T87We
model. The amplitudes and phases deviate signicantly more for
{j
0
1
} and {j
2
2
} during both seasons between these two mod-
els; again, the diurnal variation of the former and the semi-diurnal
variation of the latter are near zero at equinox in T87We. As for
CM3, it agrees well with T87We in its diurnal variation of {j
1
1
} in
December and its diminished amplitude at equinox. However, the
level of agreement is much less in all other cases where T87We
variations are not zero, particularly the amplitudes. The fact that
the agreement is best for the single local time tilted dipole coef-
cient of the list indicates that the Q
LB
ltd
Br |
and Q
|
2
s
Br |
norms.
As for the ionosphere, the CM3 parametrization has higher lati-
tudinal resolution such as to t the eld of the EEJ. Both non-local
time terms and QDconstraints allowthe CM3 eld to better respond
to the ambient eld, especially in terms of the conductivity distribu-
tion. Furthermore, the Q
Jeq
constraint injects known physical lim-
its on the conductivity patterns at night time. The GSFC(8,95-SqM)
model includes no such features. Rigid contraction and expansion
of the ionospheric eld in response to solar activity is also found in
the CM3 model, but not GSFC(8,95-SqM).
The magnetospheric eld of CM3 differs in two major ways from
that of GSFC(8,95-SqM): rst, it has the possibility of displaying
smaller-scale features in both latitude and longitude; and secondly,
it contains non-local time terms. Unfortunately, because of separa-
bility problems between the magnetospheric and ionospheric elds
inherent in the data sets used, the full impact of these differences
could not be explored.
Induced elds from time-varying external elds are included in
both models. For GSFC(8,95-SqM), the eld parametrization is
explicit, that is, independent of that of the primary sources. For CM3,
the induced eld parametrization is coupled with that of the primary
sources through an a priori conductivity model. The reason for this
is to reduce the size of the already large parameter set. Though
the independent approach is inevitably of more interest, the coupled
approach does allowfor very complicated conductivity distributions
to be included in the model, via the Q matrix of eq. (12), with no
additional parameters or computational cost.
The parametrization of elds from ionospheric coupling currents
is completely absent in the GSFC(8,95-SqM) model. This is ob-
viously a very signicant part of the measured eld for satellites
moving through these current regions, as can be attested to by Y-
component plots of Magsat dusk data (e.g. Fig. 9). In addition, the
CM3 representation is fairly sophisticated in that it includes sea-
sonal variation and a mechanism for conforming to the existing
conductivity patterns through the use of QD constraints.
7.2 Future work
While the previous section discussed in some sense how far the
comprehensive modelling effort has come, this section will hope-
fully give some clue as to how far it has to go, at least in the near
future. Much of the future work will focus on data issues: the inclu-
sion of vector data from satellite missions such as rsted, CHAMP
and others are expected to greatly enhance the validity of the models
in local times other than dawn and dusk. Furthermore, these data
will aid in the separation of ionospheric and magnetospheric elds.
One could also expand the scope of future models by considering
measurements from more magnetically disturbed times.
As for parametrization issues: the extension of the main eld
SV domain to the rsted/CHAMP epoch, the inclusion of a J