A Comprehensive Model of The Quiet-Time, Near-Earth Magnetic Field: Phase 3

Download as pdf or txt
Download as pdf or txt
You are on page 1of 37

Geophys. J. Int.

(2002) 151, 3268


A comprehensive model of the quiet-time, near-Earth
magnetic eld: phase 3
Terence J. Sabaka,
1
Nils Olsen
2
and Robert A. Langel
3
1
Raytheon ITSS, Lanham, MD 20706, USA. E-mail: [email protected]
2
Danish Space Research Institute, Juliane Maries Vej 30, DK-2100 Copenhagen , Denmark
3
Geodynamics Branch, Code 921, NASA/GSFC, Greenbelt, MD 20771, USA
Accepted 2002 February 15. Received 2001 September 28; in original form 2000 September 6
SUMMARY
The near-Earth magnetic eld is caused by sources in the Earths core, ionosphere, magneto-
sphere, lithosphere and fromcoupling currents between the ionosphere and the magnetosphere,
and between hemispheres. Traditionally, the main eld (lowdegree internal eld) and magneto-
spheric eld have been modelled simultaneously, with elds fromother sources being modelled
separately. Such a scheme, however, can introduce spurious features, especially when the spa-
tial and temporal scales of the elds overlap. A new model, designated CM3 (Comprehensive
Model: phase 3), is the third in a series of efforts to coestimate elds from all of these sources.
This model has been derived from quiet-time Magsat and POGO satellite and observatory
hourly means measurements for the period 19601985. It represents a signicant advance in
the treatment of the aforementioned eld sources over previous attempts, and includes an ac-
counting for main eld inuences on the magnetosphere, main eld and solar activity inuences
on the ionosphere, seasonal inuences on the coupling currents, a priori characterization of
the inuence of the ionosphere and the magnetosphere on Earth-induced elds, and an explicit
parametrization and estimation of the lithospheric eld. The result is a model that describes
well the 591 432 data with 16 594 parameters, implying a data-to-parameter ratio of 36, which
is larger than several popular eld models.
Key words: Earths magnetic eld, electromagnetic induction, geomagnetic variation, iono-
sphere, lithosphere, magnetosphere.
1 I NTRODUCTI ON
The terrestrial magnetic eld is composed of contributions from
many sources. Resolution of these constituent elds is of primary
importance in understanding the physical processes responsible for
their existence. This paper is concerned with the modelling of the
eld from measurements at the Earths surface and extending to
about 2000 km in altitude above that surface, a region that will
be referred to as near-Earth. Before embarking on this inverse
problem, however, a brief background into the nature of these source
phenomena is in order.
1.1 Near-Earth magnetic eld contributions
A schematic of the current sources contributing to the near-Earth
magnetic eld is given in Fig. 1. By far the most dominant of these
elds is of core origin, accounting for over 97 per cent of the eld
observed at the Earths surface and ranging in intensity from about
30 000 nT at the equator to about 50 000 nT at the poles. Accord-

Deceased , 2000 February 9.


ing to geodynamo theory, inductive interactions between the uid
motion of the liquid outer core and the geomagnetic eld not only
modify the source current so as to induce secular variation (SV)
of the eld, but sustain it against long-term decay caused by mag-
netic diffusion and Ohmic dissipation of the source current (see
Hollerbach 1996). The timescale of this eld is of the order of
centuries with typical SV magnitudes at the coremantle bound-
ary (CMB) of a few thousand nT yr
1
(see Jackson et al. 2000).
The solar-quiet (Sq) magnetic eld variation is a manifestation of
an ionospheric current system. Heating on the dayside and cooling
on the nightside of the atmosphere generates tidal winds that drive
ionospheric plasma against the geomagnetic eld, inducing electric
elds and currents in the dynamo region between 100140 km in
height. The current systemremains relatively xed to the EarthSun
line and produces regular, broad-scale daily variations that are seen
directly in the magnetograms of geomagnetic quiet days, hence
the name Sq. On disturbed days there is an additional variation,
which includes superimposed magnetic storm signatures. Because
the geomagnetic eld is strictly horizontal at the dip equator, there
is an enhancement of the effective conductivity, which results in an
enhanced eastward current, called the equatorial electrojet (EEJ),
owing along a few hundred kilometre wide band centred on the
32
C
2002 RAS
Comprehensive model: phase 3 33
Figure 1. Schematic of current sources contributing to the near-Earth magnetic eld and data sources considered in this study.
dayside dip equator (see Onwumechilli & Ozoemena 1989). In ad-
dition, auroral electrojets (AEJ) ow in the auroral belt and vary in
amplitude with different levels of magnetic activity. At the Earths
surface, the Sq elds are of the order of 1050 nT, depending upon
component, latitude, season, solar activity and time of day; the mag-
netic signature of the EEJ can be about 510 times that of Sq; and
that of the AEJ can vary widely from a few 10s nT during quiet
periods to several thousand nT during major magnetic storms (see
Volland 1984).
The eld originating in the Earths magnetosphere is a result pri-
marily of the ring-current and the currents on the magnetopause
and in the magnetotail. Currents owing on the outer boundary of
the magnetospheric cavity, known as the magnetopause, cancel the
Earths eld outside and distend the eld within the cavity. This pro-
duces an elongated tail in the antisolar direction within which sheet
currents are established in the equatorial plane. The interaction of
these currents with the radiation belts near the Earth produces a
ring-current in the dipole equatorial plane, which partially encircles
the Earth, but achieves closure via eld-aligned currents (FACs)
into and out of the ionosphere. These resulting broad-scale elds
have magnitudes of the order of 2030 nT near the Earth during
magnetically quiet periods, but can increase to several hundred nT
during disturbed times (see Olson 1984).
If displacement currents are neglected, then the current densities
associated with these external elds are solenoidal and therefore
must ow along closed circuits. Given the complex nature of the
conductivity structure in the near-Earth region, circuit closure is
sometimes achieved through currents that couple the various source
regions. At highlatitudes, the auroral ionosphere andmagnetosphere
are coupled by currents that ow along the Earths magnetic eld
lines (see Potemra 1982). The elds from these FACs have magni-
tudes that vary with the magnetic disturbance level. However, they
are always present (of the order of 30100 nT during quiet periods
and up to several thousand nT during substorms). Fields from these
currents have been detected in surface data in the Y(east) component
of the magnetic eld at low latitudes, with difculty, but are mostly
mapped using magnetometers aboard near-Earth-orbiting satellites.
There are also currents that couple the Sq currents systems in the
two hemispheres that ow, at least in part, along magnetic eld
lines. Detection of these has been reported by Olsen (1997a) using
data from the Magsat satellite. The associated magnetic elds are
generally 10 nT or less. Finally, there exists a meridional current
system that is connected to the EEJ with upward-directed currents
at the dip equator and eld-aligned downward-directed currents at
lowlatitudes (within 15

of the dip equator). Fields fromthis current


system have been detected by magnetic measurements taken on a
rocket (Musmann & Seiler 1978) and from those taken by Magsat
(Maeda et al. 1982). In the latter case, the EEJ coupling currents
resulted in elds of about 1540 nT in the Magsat data at dusk local
time.
The lithosphere is a rheological classication for that outer layer
of the Earth, which is rigid and the crust is the petrologically distinct
upper portion of this. The lithosphere contains regions for which the
temperature is below the Curie point of magnetite and other mag-
netic minerals. As a result, it can have magnetization that is either
induced by the present-day ambient eld or frozen into the rocks at
their last time of cooling below the Curie temperature, i.e. remanent
magnetization. Fields from the lithosphere are of amplitude up to
several thousand nT at the surface and at aircraft altitude, and up to
about 20 nT at the satellite altitudes considered in this paper (see
Langel &Hinze 1998). The length-scales of lithospheric elds range
from extremely broad, e.g. continentocean magnetization contrast
(Counil et al. 1991), to very local, although only features larger than
a few hundred kilometres may be resolved at satellite altitude.
The time-varying external elds induce currents in the Earths
electrically conducting interior. These currents produce a secondary
magnetic eld, the size of which depends upon the amplitude of the
inducing eld, the mantle conductivity and the period. For realistic
C
2002 RAS, GJI, 151, 3268
34 T. J. Sabaka, N. Olsen and R. A. Langel
values of mantle conductivity and induction on a global scale, the
amplitude of the induced contribution is about 40 per cent of the
inducing amplitude for periods of a few hours, about 30 per cent for
periods of a few days and less than 20 per cent for periods of more
than 6 months (see Parkinson & Hutton 1989).
1.2 Earlier modelling efforts
Historically, elds fromthe various sources have been modelled sep-
arately, or at least not all together. Under the assumption that mea-
surements are acquired in current-free regions, models of the core,
magnetospheric, ionospheric and lithospheric elds have taken the
formof gradients of Laplacian potential functions, usually in spheri-
cal coordinates. Intheir quest tomapthe magnetic eldthroughtime,
Bloxham & Jackson (1992) and later Jackson et al. (2000), with an
extension to earlier epochs, omit elds external to satellite sampling
regions and include the effects of crustal and ionospheric elds as
noise sources while constructing smooth spherical harmonic models
of the core SVat the CMBfromsatellite, observatoryrst-difference
and survey data. However, these excluded external elds, presum-
ably of magnetospheric origin, may be coestimated with internal
elds. This approach has been used by Langel & Estes (1985a,b) to
map the main eld (the low-degree spherical harmonic contribu-
tions from both the core and the lithosphere) in 1980 from Magsat
satellite data corrected for ionospheric effects. They include an ex-
ternal eld with an associated induced contribution of spherical har-
monic Y
0
1
, the time variations of which are proportional to the D
st
index. Sabaka et al. (1997) followed a similar method in modelling
long-period variation of the main and magnetospheric elds using
measurements from satellites, observatory annual-means (OAMs)
and survey data, but with annual averages of the aa index (Mayaud
1972, 1980; Rangarajan 1989) as a proxy for the D
st
index. This
coestimation of main, magnetospheric and induced elds has also
been applied to nightside rsted satellite data, culminating in the
rsted Initial Field Model (OIFM) of Olsen et al. (2000).
Spherical harmonic models of ionospheric elds have generally
been produced separately from the other eld contributions using
data from magnetic observatories and variometer stations (see e.g.
Matsushita & Maeda 1965; Malin 1973; Winch 1981; Campbell
1989; Olsen 1997b; Schmucker 1999). Langel & Estes (1985a) re-
ported detection of Sq elds in the data fromthe POGO(Polar Orbit-
ing Geophysical Observatories, comprising the OGO-2, OGO-4 and
OGO-6 satellites) satellites. Attempts to model the EEJ effects in
satellite data directly have been carried out by Langel et al. (1993),
who rst isolate the dip-latitude-dependent elds via ltering and
then t with either zonal harmonics in dipole coordinates or other
related empirical functions.
Global models of the lithospheric eld are realizable only with
satellite data, and have taken the form of various potential eld rep-
resentations. The usual approach has been to isolate the lithospheric
eld rst by removing estimates of the main, ionospheric and mag-
netospheric elds fromthe data andthentocorrelate what is believed
to be the remaining signal (see Langel & Hinze 1998, for details on
recommended procedures). Though the exact methods may deviate
from this, some examples of studies of this type can be found in
Arkani-Hamed & Strangway (1985a,b, 1986), Arkani-Hamed et al.
(1994), Cohen & Achache (1990), Counil et al. (1991), Hamoudi
et al. (1998), Ravat et al. (1995). A natural alternative is to include
higher-degree/higher-order terms in the internal eld potential ex-
pansion to capture the lithospheric contributions. Cain et al. (1984)
followed this procedure in deriving a degree/order 29 internal eld
model from Magsat data corrected for both magnetospheric and
ionospheric effects, and later, Schmitz et al. (1989) and Cain et al.
(1989a, 1990) derive even higher-degree (50) expansions using
improved ionospheric data corrections.
Currents that couple the ionosphere with the magnetosphere ow
in the ionospheric F-region at satellite altitude and therefore pro-
duce toroidal magnetic elds (see Backus 1986). Takeda & Maeda
(1983) modelled the elds caused by meridional currents as an
F-region dynamo, but perhaps the best global model of the elds
fromthese coupling currents comes fromthe work of Olsen (1997a),
who represents them as a toroidal stream function expansion within
the Magsat sampling shell under the assumption of radial currents
only.
1.3 Comprehensive approach
Each of the studies cited in the previous section enjoy varying de-
grees of success in their ability to describe the target source eld.
They all, however, suffer at some level fromthe effects of frequency
overlap between the spectra of the various source elds, both in the
spatial and temporal domains. That is, frequency range cannot be
used to absolutely distinguish between the spectra, and so bandpass
and bandstop lters alone are doomed to either keep some of the
unwanted signal or eliminate some of the signal of interest. Note
that this is different from aliasing, which is imposed by sampling
intervals and results in signals with frequencies above the Nyquist
frequency being overlapped on to those at and belowit (Kanasewich
1981).
Evidently, additional information is needed to resolve the source
contributions to near-Earth elds in a realistic manner. One possibil-
ity is to consider the radial positions of the various source regions
relative to the available data, which is shown in the schematic in
Fig. 1. Core, lithospheric and induced elds would be internal to
both surface and satellite data while the magnetospheric eld would
be external. The ionospheric eld would be external to surface data,
but internal to satellite data. Thus, surface data could separate iono-
spheric and magnetospheric from core, lithospheric and induced
elds and satellite data could separate magnetospheric from core,
lithospheric, induced and ionospheric elds. This suggests that a
joint analysis of both surface and satellite data could theoretically
resolve parametrizations between magnetospheric, ionospheric and
internal sources, but only if the parameter set is treated consis-
tently between data types, which implies that they be coestimated.
Of course, additional information would be needed to separate the
various internal elds, e.g. space and timescale differences, a priori
mantle conductivity models, etc. This simultaneous inversion for
parameters describing all sources will be termed the comprehen-
sive approach, and models of this kind could provide the reference
elds needed in more rened studies where the separation/removal
of the various source contributions is an issue.
This paper reports on the third in a series of efforts to derive
progressively more sophisticated models using the comprehensive
approach, and many of the details omitted here are included in the
associated NASA Technical Memorandum of Sabaka et al. (2000).
The rst phase, reported by Sabaka & Baldwin (1993), culminated
in a model known as GSFC(12,93), while the second phase, re-
ported by Langel et al. (1996), culminated in the GSFC(8,95-SqM)
model. These models are based upon magnetically quiet data from
the POGOand Magsat spacecraft and OAMs and observatory hourly
means (OHMs). Because of their limited scope, the Magsat vector
data poleward of 50

dipole latitude are not used in order to avoid


the auroral FACs, while the Magsat dusk data are corrected for the
effects of the EEJ and associated meridional currents. Details of the
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 35
model parametrizations are not given here, but qualitatively, both
include representations of the main eld, its SV and the local time
(Sun-synchronous) modes of the magnetospheric and ionospheric
elds; both include ring-current variations through the D
st
index;
and both include explicit parametrizations for elds induced by the
time-varying external elds. The GSFC(8,95-SqM) also includes
seasonal variations in the magnetospheric and ionospheric elds.
The scope of the model presented in this paper is much wider
and its attention to detail much higher than its predecessors. Be-
ing the culmination of the third phase of work along these lines,
this model will be designated as CM3 (denoting Comprehensive
Model: phase 3). Its aimis to describe ner systematic detail present
in the data sets with a physically plausible model of least complex-
ity, which implies a more proper distribution of the signal amongst
previously considered sources and amongst newly parametrized
sources, i.e. coupling currents and small-scale lithosphere. In ad-
dition to these new sources, other improvements have been made
such as a smooth main eld SV; a high-resolution, high-efciency
ionospheric parametrization with physical conductivity constraints
and main eld interaction; a magnetospheric parametrization with
maineldinteraction; andanimplicit parametrizationof the induced
elds associated with the time-varying external elds via a radially
varying a priori mantle conductivity model. In the remainder of
the paper, a description of the data, parametrization and method of
estimation leading to the CM3 model will be given followed by a
discussion of the model properties in terms of consistency and phys-
ical plausibility. In conclusion, model availability will be discussed
and future directions will be outlined.
1.4 Coordinate systems
It is convenient at this point to introduce the various spatial and
temporal coordinate systems used throughout this paper. There are
basically three spatial systems employed for different parts of the
model. The usual geographic spherical polar coordinates, (r. . ),
are dened with the z-axis through the north geographic pole, x
through Greenwich and y completing the right-handed system, and
are used for internal eld representations. The dipole spherical polar
coordinates, (r.
d
.
d
), are dened with respect to a tilted dipole-
moment vector, which in this study is supplied by the GSFC(12,83)
model of Langel & Estes (1985b). The z-axis points along the mo-
ment vector out of the northern geographic hemisphere, the x-axis
is perpendicular to z in the half-plane containing the south geo-
graphic pole and y completes the right-handed system. This system
is used primarily to represent distant external (i.e. magnetospheric)
eld contributions. Finally, the quasi-dipole (QD) coordinates pro-
posed by Richmond (1995), (h
q
.
q
.
q
), are also tied to a magnetic
eld model, which in this study is the DGRF 1980 model (IAGA
Division I Working Group 1 1981). The QD longitude,
q
, is the
dipole longitude of the apex (highest point above the Earths sur-
face, given by h
A
) of the magnetic eld line passing through the
point. The QD latitude,
q
, is given by

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

) dipole latitude vector data in the analysis,


particularly the X and Y components, which are sensitive to high-
latitude FACs. Accordingly, X and Y vector components were added
at samplingpoints polewardof 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

, and therefore one may be able to


t this feature with far fewer QD symmetric functions.
Before proceeding, certain properties of the QDcoordinates must
be articulated. First, the Laplacian operator does not separate in QD
coordinates, rendering closed-form solutions out of the question.
Secondly, though the QD colatitude,
q
, and longitude,
q
, do chart
the sphere, their coordinate lines change slightly with radius, r.
With this in mind, consider a set of basis functions, which possess
the QD symmetry in two dimensions rather than three, perhaps on a
constant-coordinate surface. If a sphere is chosen, then it is natural
to think of the QD angular coordinates as the formal arguments to
the usual surface harmonics. These functions may be expanded in
terms of the surface harmonics in dipole coordinates as
P
l
k
(cos
q
) exp il
q
=
Nmax

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

. This suggests functions of the form


T
l
k.j
(r) =

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

damps the higher-degree


harmonics much more severely than Q
Jeq
, as intended.
4.1.4 Minimal magnetospheric local-time dipole deviations
It is anticipated that the magnetospheric eld expansion of eq. (29)
includes many more coefcients than can be reliably estimated from
the data at hand, especially those describing deviations froma dipole
in MLT. Experience from earlier phases of modelling suggests that
excessive cross-talk or correlations between the ionospheric and
non-D
st
-dependent magnetospheric expansions will probably exist
as a result of poor eld separation caused by limited satellite vector
data coverage in local time. Therefore, a magnetospheric solution is
sought which is smooth in some sense that will reduce this coupling.
Specically, dene a fth quadratic norm, Q
LB
ltd

, 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

centred on the geographic poles (Langel & Estes 1985a). This,


combined with the fact that J
r
is expressed in dipole coordinates,
makes damping the polar regions a necessity for both dawn and
dusk. Furthermore, since J
r
shows little structure at low and mid-
latitudes during dawn, there is no need to introduce a
d
inuence
function as in Q

2
s
J
eq. p>0

. Although the EEJ coupling currents are


present at dusk along the dip equator, the inclusion of an inuence
function that is small only at low dipole latitudes is complicated,
and so it is omitted in this study.
To this end, dene a sixth type of quadratic norm, Q
| Jr |
, which
measures the mean-square magnitude of J
r
on a sphere at Magsat
altitude (r = a +h
m
with h
m
= 450 km), O, over UT, given by
Q
| Jr |
(x) =
_
T
0
_
O
[J
r
(t.
d
.
d
)]
2
dOdt
__
T
0
_
O
dOdt
_
1
. (58)
where
J
r
(t. r) =
_
_
_
1
j
0
(a +h
m
)
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
_

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

) OHM-MUL processes are considered to be


isotropic; the high dipole latitude (|90


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

column of Table 3 (see Section 5).


This procedure also provides an opportunity to reject gross out-
liers with respect to the A model. Specically, residuals in the
Magsat mid/low and high dipole latitude data sets greater that 25
and 100 nT, respectively, are rejected as are OHM-MUL residuals
greater than 150 nT and POGO residuals greater than 25 nT. The
resulting measurement counts for the various data sets are listed in
the number column of Table 3 (see Section 5).
4.3 Application
Two iterations of the Gauss method were required from an x
0
pro-
vided by a preliminary model known as B, a successor to A, which
is described more fully by Sabaka et al. (2000). This B model,
also referred to as the GSFC/CU(12,96) model, has been used by
Purucker et al. (1997) to study northsouth-trending anomalies of
lithospheric origin, particularly over Australia.
In the present study, the number of observations (dimd) is
591 432, the number of model parameters (dimx) is 16 594 and
the number of quadratic smoothing norms is seven (x
a
i
= 0 for
each).
Because seven norms are employed in this study and because
the number of parameters is large, a typical trade-off study be-
tween the various terms in the least-squares objective function is
impractical and so the more qualitative method of visual inspection
is used. This is thought to sufce for two reasons: rst, the norms
measure physically meaningful quantities about which qualitative
knowledge abounds and which can be checked visually; and sec-
ondly,
i
have a logarithmic inuence on the norms and so only
an order of magnitude accuracy is needed. The additional criterion
of low correlation between ionospheric and magnetospheric coef-
cients was instrumental in selecting the Q
LB
ltd

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

W at noon MLT. This panel shows a


progression of residuals from top to bottom, 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
. Note the clear EEJ signature in the residuals of the second
member, which the model is able to reproduce. Though the t is
satisfactory for most of the pass (8 nT), it begins to steadily di-
verge north of about 20

geographic latitude, reaching a maximum


magnitude of 7.7 nT at 51

. This may be a result of the inuence


of external current systems, particularly unmodelled ring-current
dynamics.
The t to the Z component for Magsat dusk pass 263, which
was not included in the model, is shown in a similar format in
the bottom panel of Fig. 8. Again, the total t is satisfactory for
practically all of the pass (15 nT). The corresponding X- and Y-
component residual progressions are shown in Fig. 9. The ts are
within about 15 nT for all of X and most of Y. The obvious ex-
ception is found in the polar region south of about 60

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 at noon while the Magsat pass crosses at 129

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

; geometrically, the mid-point between being fully cor-


related and uncorrelated.
Three of the ve major categories of correlations at this level
involve spatial and temporal separability problems within the same
eld source and indicate overparametrizations: in the main eld
SV, where B-spline density per harmonic is high; in the iono-
spheric E-region, where the latitudinal resolution is high; and in the
ionospheric F-region coupling currents, particularly the seasonal
variation.
A fourth category includes correlations amongst OHM biases, at
the same and different locations. Correlations across breaks in a par-
ticular station record may be caused by both segments cross-talking
with the main eld SV. Correlations between biases at spatially dis-
tinct points are intriguing because they may provide insight into
spatial correlation lengths of crustal eld sources on local scales.
To investigate this, a global map was produced in which the loca-
tions of spatially distinct observatories were connected by a line
if any bias component of one was correlated with any bias com-
ponent of the other, with an absolute value above 0.7 (see the top
panel of Fig. 11). The box outlines the European sector (exploded
in the bottom panel) and the circles have been added to help guide
the eye to some of the more obscure lineations. First, note that none
of the line segments is longer than about 6

of arc and most are


much shorter. This corresponds well with the resolution limit of the
N
max
= 65 truncation level of the core and lithospheric eld ex-
pansion, and one would expect that any coherence at lengths longer
than this would be described by the model. Secondly, note that the
largest concentration and length of lines are found in Europe. This
follows from the fact that the highest concentration of observato-
ries is found in Europe. If a signicant crustal correlation length is
present elsewhere in the world, but is only sampled by one obser-
vatory, then this exercise will fail to detect it. Although a detailed
interpretation will not be attempted here, it must be mentioned that
some of the lineation patterns in Europe do appear to be related
to known geological features. The polygon dened by the Troms
(TRO), Abisko (ABK) and Kiruna (KIR) observatories agrees well
with the location of the Kiruna crustal anomaly, while the eastwest
lineation between Nurmij arvi (NUR) and Voyeykovo (LNN) may
be related to the Kursk anomaly. The complex of F urstenfeldbruck
(FUR), Wien Kobenzl (WIK), O Gyalla Pesth (OGY) and Tihany
(THY) seems to follow the fabric of the Alpine region, while the
eastwest lineations in the Iberian peninsula may also be following
regional trends.
The last category represents the only signicant cross-talk be-
tween different eld sources. It involves negative correlations be-
tween the parameters of the non-D
st
-dependent magnetosphere and
the ionospheric E-region with similar spatial and temporal be-
haviour. In a preliminary CM3-type model that did not include the
Q
LB
ltd

smoothing, these correlations exist between harmonics hav-


ing broad-scale latitudinal variation. For CM3, these correlations ex-
ceed the 0.7 threshold only for the tilt component (k = l = p = 1)
of the MLTdipole, andmost of that is inthe noonmidnight direction
(real part).
While a direct inspection of R
x
is certainly a worthwhile exer-
cise, it can be misleading: most physically meaningful quantities
are functions of many parameters, and even though correlations
between any given pair may be negligible, the cumulative affect
may be signicant. It is natural then to wonder whether any further
cross-talk between eld sources does exist, and so global correlation
C
2002 RAS, GJI, 151, 3268
54 T. J. Sabaka, N. Olsen and R. A. Langel
|Observatory Bias Correlations| > 0.7
180 210 240 270 300 330 0 30 60 90 120 150 180
-90
-60
-30
0
30
60
90
0 30
30
60
ABK
KIR
TRO
ALM
SFS
BEL
SWI
BFE
NGK
RSV
COI
TOL
EBR
ESK
STO
FUR
WIK
HAD
VAL
MMK
LOP
NUR LNN
OGY
THY
STP
ODE
WNG
WIT
Figure 11. Observatory bias correlations between spatially distinct locations with absolute values above 0.7 shown as connecting lines (cylindrical equidistant
projection). This includes correlations between any component of one observatory bias with any component of another. Circles have been added to the top
panel as a visual aid in locating some of the more obscure lineations.
maps between like pairs of vector eld components from suspected
sources were generated. These include magnetosphere with iono-
sphere (Fig. 12) and lithosphere with ionosphere (Fig. 13). Correla-
tions between elds from coupling-currents and potential elds and
correlations between lithospheric and magnetospheric elds were
found to be negligible. Since it is of interest to gauge the worst-
case separability problems, these maps are generated at the summer
solstice when Magsat data, probably the most dominant separating
agent for the pairs of interest, are absent. Of the two pairs of eld
sources, the correlations betweenmagnetosphere andionosphere are
stronger, yet the largest correlation magnitude is only 0.41, which
translates to 66

. Therefore, eld separability in CM3 is judged


to be satisfactory for most applications. Fig. 12 shows ranges of
(0.16. 0.21), (0.02. 0.41), (0.29. 0.09) for XX, YY and Z Z,
respectively. Since the maps are centred on noon MLT, one can see
that the most dominant features appear to be Earth-xed. This may
indicate a slight cross-talk between the magnetosphere and the static
portion of the ionosphere introduced as a baseline for nightside cur-
rent minimization. It is also interesting to notice that on average YY
is strongly positive, Z Z is strongly negative and XX is mixed.
For the lithospheric and ionospheric correlations, Fig. 13 shows
ranges of (0.16. 0.06), (0.26. 0.06), (0.29. 0.04) for XX, YY
and Z Z, respectively, which are predominantly negative. Note that
the colour scale has nowchanged fromFig. 12 in order to accentuate
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 55
Figure 12. Correlation maps of like vector components from magnetospheric and ionospheric elds at r = 6821.2 km on 1980 June 21 centred on noon MLT,
but for different MUTs (Mollweide projection).
Figure 13. Correlation maps of like vector components from lithospheric and ionospheric elds at r = 6771.2 km on 1980 June 21 centred on noon MLT,
but for different MUTs (Mollweide projection).
C
2002 RAS, GJI, 151, 3268
56 T. J. Sabaka, N. Olsen and R. A. Langel
the weaker structure. Most features appear to be more closely xed
in MLT than in the previous case, and the strongest features are
in the northern polar regions of XX and Z Z. The more intense
localized spots in the Pacic basin of the XX and Z Z maps do not
appear to correlate with actual features in the lithospheric anomaly
maps. However, it must be reiterated that the correlations between
the magnetosphere and the ionosphere, and between the lithosphere
and the ionosphere are weak at best.
5.3 Resolution and calibration
Several quantities related to parameter resolution and error calibra-
tion have been computed and are listed in Table 4. They reect both
the data subsets dened in Table 3 and the various norms discussed
in Section 4.1. These quantities include two additional classes of
mists: data subset mists, s
2
d
i
, which measure the efciency of the
model in tting data in that subset; and prior mists, s
2
a
i
, which mea-
sure the departure of the model from its ith a priori preferred state
per DOF in that norm. They are dened as
s
2
d
i
=
[d
i
a
i
( x)]
T
W(d
i
a
i
( x))
N
i
tr
_
R
d
i
_ . (62)
s
2
a
i
=

i
Q
i
( x)
M
i
tr
_
R
a
i
_. (63)
with
R
d
i
= C
x
A
T
i
WA
i
. (64)
R
a
i
=
i
C
x
A
i
. (65)
where N
i
is the number of observations in the ith data subset,
M
i
= dimx
a
i
, and A
i
is the submatrix of A associated with the
Table 4. CM3 resolution and calibration information. L
r
i
and L
e
i
are the weighted residual and
error variances, respectively; R
d
i
and R
a
i
are the data and norm resolution matrices, respectively;
N
i
and M
i
are the dimensions of the data and preferred model state vectors, respectively; s
2
d
i
and s
2
a
i
are the data and prior mists, respectively; and tr[] and E[] are the trace and expectation operators,
respectively.
Data subset L
r
i
( x) N
i
tr[R
d
i
] E[L
r
i
( x)] s
2
d
i
OHM-1AM 94 738.0 81 712 924.9 80 787.1 1.17
OHM-MUL
|90

d
| - 50

115 849.2 166 334 1527.6 164 806.4 0.70


|90

d
| 50

194 267.7 197 422 1060.5 196 361.5 0.99


Magsat dusk 23 663.1 39 488 2672.4 36 815.6 0.64
Magsat polar dusk 12 251.0 15 973 212.0 15 761.0 0.78
Magsat dawn 27 327.1 44 136 2572.1 41 563.9 0.66
Magsat polar dawn 14 861.5 16 928 191.9 16 736.1 0.89
POGO decimated 17 367.4 22 685 461.2 22 223.8 0.78
POGO pass 5018.6 6 754 133.5 6,620.5 0.76
Subtotal 505 343.6 591 432 9756.2 581 675.8 0.87
Norm L
e
i
( x) M
i
tr[R
a
i
] E[L
e
i
( x)] s
2
a
i
Q
|

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

1149.3 4910 1448.9 3461.1 0.33


Q
LB
ltd

2907.6 740 612.9 127.1 22.88


Q
| Jr |
dawn 743.8 1044 474.4 569.6 1.31
Q
| Jr |
dusk 744.6 1044 524.2 519.8 1.43
Subtotal 49 288.0 18 328 6837.8 11 490.2 4.29
Grandtotal 554 631.5 609 760 16 594.0 593 166.0 0.94
ith data subset. In the table, L
r
i
( x) and L
e
i
( x) refer to the numera-
tors of eqs (62) and (63), respectively, and E[] is the expectation
operator. From the fourth column, which gives the trace of the res-
olution matrix corresponding to either the ith data subset, R
d
i
, or
the norm, R
a
i
, it can be seen that the data are resolving about 59
per cent of the 16 594 total parameters. Of these, about 36, 58 and
6 per cent are resolved by the observatories, Magsat, and POGO,
respectively. As for the norms, about 13, 0.3, 10, 9 and 4 per cent of
the parameters are resolved by Q
|

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

norms, suggesting that


their inuence be signicantly decreased. However, since L
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

, at which the spectra become level may


be an indication of the maximum depth of the current source layers
associated with that part of the magnetic eld, and for a power law
of the form R
n
= c(d)
n
at radius R it can be shown that (Cain et al.
1989b)
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

E). Other CM3 north


south lineations, such as those in the South American mid-continent
region just south of the dip equator, may be manifestations of exter-
nal eld contamination, though it is unclear at this point.
To facilitate the comparison at high latitudes, B
r
polar maps above
60

N and below 60

S are shown in Fig. 19. Because the RLPAA


map is derived from Magsat data alone, it is not considered valid
above about 83

N or below about 83

S. However, within the region


of Magsat coverage, the maps agree very well in general shape and
location of the major known anomalies. Again, the anomaly inten-
sities run generally higher for the CM3 maps. There are a couple of
notable differences, such as the slight shift of the northern Greenland
anomaly (I) from a position on the northwest coast in RLPAA to a
more north-central position in CM3, and the emergence of a sig-
nicant negative anomaly centred over the south magnetic dipole
position (J) (79

S, 109

E) in CM3 that is absent in RLPAA. The


C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 61
Figure 18. Global map of the CM3 lithospheric contribution (degrees 1542) to B
r
on the sphere r = 6771.2 km (cylindrical equidistant projection). Letters
identify particular anomalies referenced in the discussion.
Figure 19. Polar maps above 60

N (top) and below 60

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

, respectively. Thus, the gap


in polar coverage for CM3 has been reduced to caps of about 3

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

towards MLT noon, i.e. {j


1
1.0.1
} =1.0 nT, the seasonal
uctuation about this mean is much larger, being predominantly in
the annual term, with the peak tilt magnitude being 15.3

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

smoothing, denoted CNXS; both for


magnetospheric terms alone and for the sum of magnetospheric and
ionospheric terms of similar indices. A value of F
10.7
= 140.0
10
22
Wm
2
Hz
1
was usedtocompute the ionospheric coefcients,
and all coefcients are in geographic spherical-polar coordinates.
Comparing coefcients of the form {j
m
np.g
+c
m
np.g
}, it can be seen
that the phases, , are in fair agreement (within about 40

) 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

smoothing, which is damping


Table 7. Magnetospheric annual variation (nT) from Malin & Mete Isikara
(1976) compared with CM3-type models. A value of F
10.7
= 140.0
10
22
W m
2
Hz
1
was used to compute the ionospheric coefcients, and
all coefcients are in geographic spherical-polar coordinates.
{c
m
np
(t )} M&MI CM3 CNXS
A A A
{j
0
1.0.g
+c
0
1.0.g
} 3.8 190

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

Equinox 1.4 357

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

Equinox 0 0 0.1 338

2.6 357

{j
1
1
} December 0.8 181

1.1 181

0.1 171

0.7 177

Equinox 0.7 174

0.7 174

0.1 230

0.6 151

much of the amplitude of the n = 2 terms in the CM3 model. This


is clearly seen by comparing coefcients of the form {j
m
np.g
}; the
amplitudes of {j
0
2.0.g
} and {j
1
2.1.g
} are an order of magnitude
less in CM3. If M&MI are correct, then perhaps this damping is
excessive, which agrees with the conclusions of Section 5.3.
Olsen (1996) has investigated possible magnetospheric eld con-
tributions to daily magnetic eld variations. He examined several
semi-empirical models of the magnetospheric eld, reporting pri-
marily on results from Tsyganenko (1987, 1989) and, in particular,
the model designated T87We in Peredo et al. (1993). If the real part
of a particular coefcient, say {j
m
n
(t
s
. t
mut
)}, is Fourier analysed
for diurnal ( p =1) or semi-diurnal ( p =2) periodicities in MUT at
a xed season (specied by t
s
) such that

_
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

smoothing may again


be excessive. Hence, the better agreement of the diurnal varia-
tion of {j
0
1
} and the semi-diurnal variation of {j
2
2
} at equinox
for CM3 compared with the zero values for T87We is considered
fortuitous.
The magnetospheric model predictions of the scalar intensity
along the POGO pass shown in Fig. 8 and the vector components
along the Magsat dusk pass 263 in Figs 8 and 9 are clearly those
of an axial dipole whose moment is aligned, on average, nearly per-
pendicular to the orbit normals of the satellites. Furthermore, the
behaviour of the model as a function of the D
st
index is illustrated in
the X and Z components of Tucson in 1967 (see Fig. 5) where very
large discontinuities are found in the magnetic records between the
quietest days chosen for May and June. They are evidently the result
of a ring-current adjustment, which is described quite well by D
st
,
and hence, the model ts them with little problem.
6.6 Fields from ionospheric coupling currents
The radial current density (J
r
) of the F-region coupling currents
from CM3 at dusk local time is shown in Fig. 12 for March and
December 21. Note that this gure does not show J
r
at any par-
ticular UT as a function of longitude, but rather, it shows the J
r
predicted to ow at different longitudes at dusk local time. As ex-
pected, the largest radial current ows at polar latitudes. During the
morning (not shown), currents ow into the ionosphere (J
r
- 0)
at the poleward boundary of the polar oval (region 1 currents), and
out of the ionosphere at the equatorial boundary (region 2 currents).
During the evening the current direction is reversed. The evening
data also show upward currents at the dip equator and downward
currents at nearby low latitudes. This is the radial component of
the meridional current system of the EEJ, rst observed in Magsat
data by Maeda et al. (1982). There is no evidence for such a current
system in the morning sector. A strong seasonal variation can be
seen in J
r
between March and December. During the southern sum-
mer it is much more intense in the southern polar oval region than
in the northern one, while during equinox the intensities are more
equal. In addition, the radial component of the meridional current
system is better dened during equinox than during solstice. Al-
though these ndings are in keeping with those of Olsen (1997a),
he nds current intensities that are in general signicantly higher
than those predicted by CM3. This discrepancy, however, may be
attributed to the difference in data selection: all days versus quiet
days.
Olsen (1997a) has also detected larger-scale upward currents in
the north and downward currents in the south in the evening, and
opposite this in the morning, for the December Magsat data. These
are obscure in Fig. 22, if present at all, and may probably be a
result of the application of the Q
| Jr |
smoother at both the dawn and
dusk local times, and again to data selection. Such interhemispheric
coupling currents are small or absent during the equinoxes and are
expected to reverse during the northern summer.
The toroidal eld signature of the meridional current system may
be clearly seen in the Y component of Magsat dusk pass 263 in
Fig. 9. Because this pass occurs on 1979 November 19, near the
southern summer, the current density vortex just to the south of
the dip equator is stronger than that just to the north, resulting in a
stronger eastward and weaker westward eld just south and north
of the dip equator, respectively. The CM3 model, in turn, predicts
this asymmetry very well. The high-amplitude excursions in the X
and particularly the Y components at high latitudes are probably also
of the eld-aligned variety. However, these are also probably very
transient in nature, and are thus tted very poorly by a model of this
type, which sees mean temporal and spatial eld effects best.
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 65
Figure 22. Global maps of the radial component of the ionospheric coupling currents (J
r
in eq. 59) at dusk on the sphere r = 6821.2 km for March 21 and
December 21 (cylindrical equidistant projection).
C
2002 RAS, GJI, 151, 3268
66 T. J. Sabaka, N. Olsen and R. A. Langel
7 CONCLUS I ONS
The paradigm of comprehensive modelling is quite worthwhile, yet
quite formidable. Progress has been slow and incremental, and with
new satellite missions and new technology on the horizon, a nal
model is simply an ideal to work towards. In this nal section an
attempt will be made to gauge the position of CM3 and its method-
ology in the continuum of comprehensive modelling.
7.1 New features
Perhaps the best way to assess the progress made in the current phase
is to simply compare it with the previous phase. Hence, features of
the CM3 model which are newwith respect to the GSFC(8,95-SqM)
model (Langel et al. 1996) will be enumerated in this section.
The main eld and attendant SV portions (triple summation in
eq. 5) of the models are identical, but CM3 includes a static repre-
sentation of the high-degree lithospheric eld (double summation
in eq. 5), which has successfully captured most of the known crustal
anomaly elds seen at satellite altitude. In addition, the SV is now
much more plausible because of the regularization imparted by the
Q
|

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

component to the coupling current density, and a continuous diur-


nally varying toroidal eld will be pursued. More long-term goals
include the proper treatment of elds from the ionosphere and as-
sociated coupling currents with respect to polar currents that are
dependent upon the interplanetary magnetic eld, the extension to
more magnetically disturbed periods, and the incorporation of more
complicated a priori conductivity models.
7.3 Possible uses
The success of comprehensive modelling is in part driven by its util-
ity to the scientic community. The method of coestimating elds
from several sources and its affect on model consistency is of sci-
entic interest in its own right; however, additional merit of the
comprehensive models lies in their use as application tools, or refer-
ence models. Indeed, with the possible exception of the high-degree
lithospheric eld where new, physically meaningful features might
be found, most source elds are parametrizedsoas tomodel the well-
known, regular, quiet-time features. Hence, comprehensive models
are well qualied for removing known elds from the data so as to
not obfuscate that which is unknown. However, one must be cautious
when applying them outside their scope, that is, extrapolating them
to regimes that were not sampled by the data sets used in deriving
them. Uncritical application of CM3 during magnetically disturbed
conditions, for instance, is dangerous.
7.4 Availability
In accordance with the previous section, a forward modelling code is
available, which predicts the various CM3 source elds given spatial
and temporal positions and D
st
and F
10.7
values. This code is in the
form of an ANSI standard Fortran subroutine called CM3FIELD.
It returns the local (north, east, down) components of B in nT on a
sphere or on the IAU1966 ellipsoid from internal, primary and in-
duced ionospheric, primary and induced magnetospheric, and dawn
or dusk coupling current eld sources. Two evaluations of the in-
ternal eld are accommodated per two given spherical harmonic
degree ranges. This is helpful when separate predictions are desired
from the core-dominated and lithosphere-dominated portions of the
expansion. The CM3 model information is input to the CM3FIELD
subroutine on the initial call from a sequentially accessed ASCII
le. Both the forward code le and the CM3 model information le
are available from the authors by request.
ACKNOWLEDGMENTS
We would like to thank Art Richmond for his expertise in several
external-eld-related issues and for help with the initial stages of
the manuscript. We also thank Richard Holme and Gauthier Hulot
for their thorough reviews. All gures were produced with GMT
C
2002 RAS, GJI, 151, 3268
Comprehensive model: phase 3 67
(Wessel & Smith 1991). TJS was supported under NASA Contract
NAS5-31760.
REFERENCES
Arkani-Hamed, J. & Strangway, D.W., 1985a. Intermediate-scale magnetic
anomalies of the Earth, Geophysics, 50, 28172830.
Arkani-Hamed, J. & Strangway, D.W., 1985b. Lateral variations of apparent
susceptibility of lithosphere deduced from Magsat data, J. geophys. Res.,
90, 26552664.
Arkani-Hamed, J. &Strangway, D.W., 1986. Band-limited global scalar mag-
netic anomaly map of the Earth derived fromMagsat data, J. geophys. Res.,
91, 81938203.
Arkani-Hamed, J., Langel, R.A. &Purucker, M.E., 1994. Magnetic anomaly
maps of Earth derived from POGO and Magsat data, J. geophys. Res., 99,
24 07524 090.
Backus, G.E., 1986. Poloidal and toroidal elds in geomagnetic eld mod-
eling, Rev. Geophys., 24, 75109.
Barton, C.E., 1997. International geomagnetic reference eld: the seventh
generation, J. Geomag. Geoelectr., 49, 123148.
Bloxham, J. & Jackson, A., 1992. Time-dependent mapping of the magnetic
eld at the coremantle boundary, J. geophys. Res., 97, 19 53719 563.
Bloxham, J., Gubbins, D. & Jackson, A., 1989. Geomagnetic secular varia-
tion, Phil. Trans. R. Soc. Lond., A., 329, 415502.
Cain, J.C., Schmitz, D.R. & Muth, L., 1984. Small-scale features in the
Earths magnetic eld observed by Magsat, J. geophys. Res., 89, 1070
1076.
Cain, J.C., Wang, Z., Kluth, C. & Schmitz, D.R., 1989a. Derivation of a
geomagnetic model to n = 63, Geophys. J., 97, 431441.
Cain, J.C., Wang, Z., Schmitz, D.R. & Meyer, J., 1989b. The geomagnetic
spectrum for 1980 and corecrustal separation, Geophys. J., 97, 443447.
Cain, J.C., Holter, B. & Sandee, D., 1990. Numerical experiments in geo-
magnetic modeling, J. Geomag. Geoelectr., 42, 973987.
Campbell, W.H., 1989. The regular geomagnetic eld variations during quiet
solar conditions, in Geomagnetism, Vol. 3, pp. 385460, ed. Jacobs, J.A.,
Academic Press, London.
Chapman, S. & Bartels, J., 1940. Geomagnetism,Clarendon Press, Oxford.
Cohen, Y. & Achache, J., 1990. New global vector magnetic anomaly maps
derived from Magsat data, J. geophys. Res., 95, 10 78310 800.
Counil, J., Cohen, Y. & Achache, J., 1991. The global continentocean mag-
netization contrast, Earth planet. Sci. Lett., 103, 354364.
Hamoudi, M., Cohen, Y. & Achache, J., 1998. Can the thermal thickness
of the continental lithosphere be estimated from Magsat data?, Tectono-
physics, 284, 1929.
Hollerbach, R., 1996. On the theory of the geodynamo, Phys. Earth planet.
Inter., 98, 163185.
Holme, R. & Bloxham, J., 1996. The treatment of attitude errors in satellite
geomagnetic data, Phys. Earth planet. Inter., 98, 221233.
IAGA Division I Working Group 1, 1981. International Geomagnetic Ref-
erence Fields: DGRF 1965, DGRF 1970, DGRF 1975, and IGRF 1980,
EOS, Trans. Am. geophys. Un., 57, 120121.
Jackson, A., 1994. Statistical treatment of crustal magnetization, Geophys.
J. Int., 119, 991998.
Jackson, A., Jonkers, A.R.T. & Walker, M.R., 2000. Four centuries of ge-
omagnetic secular variation from historical records, Phil. Trans. R. Soc.
Lond., A., 358, 957990.
Kaiser, J.F., 1974. Nonrecursive digital lter design using the Io-sinh window
function, Proc. IEEE Int. Symp. Circuits Syst., 1974, 2023.
Kanasewich, E.R., 1981. Time sequence analysis in geophysics, pp. 1480,
The University of Alberta Press, Edmonton.
Langel, R.A., 1987. The main geomagnetic eld, in Geomagnetism, Vol. 1,
pp. 249512, ed. Jacobs, J.A., Academic Press, London.
Langel, R.A., 1990. Global magnetic anomaly maps derived from POGO
spacecraft data, Phys. Earth planet. Inter., 62, 208230.
Langel, R.A. & Baldwin, R.T., 1991. Geodynamics branch data base for
main magnetic eld analysis, NASA Tech. Mem., 104542, 1101.
Langel, R.A. & Estes, R.H., 1982. A geomagnetic eld spectrum, Geophys.
Res. Lett., 9, 250253.
Langel, R.A. & Estes, R.H., 1985a. Large-scale, near-Earth magnetic elds
from external sources and the corresponding induced internal eld,
J. geophys. Res., 90, 24872494.
Langel, R.A. & Estes, R.H., 1985b. The near-Earth magnetic eld at 1980
determined from Magsat data, J. geophys. Res., 90, 24952510.
Langel, R.A. & Hinze, W.J., 1998. The Magnetic Field of the Earths Litho-
sphere: The Satellite Perspective, pp. 1429, Cambridge University Press,
Cambridge.
Langel, R.A., Coles, R.L. & Mayhew, M.A., 1980. Comparisons of mag-
netic anomalies of lithospheric origin measured by satellite and airborne
magnetometers over western Canada, Can. J. Earth. Sci., 17, 7.
Langel, R.A., Estes, R.H. & Mead, G.D., 1982. Some new methods in ge-
omagnetic eld modeling applied to the 19601980 epoch, J. Geomag.
Geoelectr., 34, 327349.
Langel, R.A., Purucker, M.E. &Rajaram, M., 1993. The equatorial electrojet
and associated currents as seen in Magsat data, J. Atmos. Terr. Phys., 55,
12331269.
Langel, R.A., Sabaka, T.J., Baldwin, R.T. & Conrad, J.A., 1996. The near-
Earth magnetic eld from magnetospheric and quiet-day ionospheric
sources and how it is modeled, Phys. Earth planet. Inter., 98, 235
267.
Lowes, F.J., 1974. Spatial power spectrum of the main geomagnetic eld,
and extrapolation to the core, Geophys. J. R. astr. Soc., 36, 717730.
Maeda, H., Iyemori, T., Araki, T. & Kamei, T., 1982. New evidence of a
meridional current system in the equatorial ionosphere, Geophys. Res.
Lett., 9, 337340.
Malin, S.R.C., 1973. Worldwide distribution of geomagnetic tides, Phil.
Trans. Roy. Soc. Lond., A., 274, 551594.
Malin, S.R.C. & Gupta, J.C., 1977. The Sq current system during the Inter-
national Geophysical Year, Geophys. J. R. astr. Soc., 49, 515529.
Malin, S.R.C. &Mete Isikara, A., 1976. Annual variation of the geomagnetic
eld, Geophys. J. R. astr. Soc., 47, 445457.
Matsushita, S. & Maeda, H., 1965. On the geomagnetic quiet solar daily
variation eld during the IGY, J. geophys. Res., 70, 25352558.
Mauersberger, P., 1956. Das Mittel der Energiedichte des geomagnetischen
Hauptfeldes an der Erdober ache und seine s akulare

Anderung, Gerlands
Beitr. Geophys., 65, 207215.
Mayaud, P.N., 1972. The aa indices: a 100-year series characterizing mag-
netic activity, J. geophys. Res., 77, 6870.
Mayaud, P.N., 1980. Derivation, meaning and use of geomagnetic indices,
Geophys. Monogr. Am. geophys. Un., 22.
Musmann, G. & Seiler, E., 1978. Detection of meridional currents in the
equatorial ionosphere, Geophys. J., 44, 357372.
Olsen, N., 1993. The solar cycle variability of lunar and solar daily geomag-
netic variations, Ann. Geophys., 11, 254262.
Olsen, N., 1996. Magnetospheric contributions to geomagnetic daily varia-
tions, Ann. Geophys., 14, 538544.
Olsen, N., 1997a. Ionospheric F-region currents at middle and low latitudes
estimated from Magsat data, J. geophys. Res., 102, 45634576.
Olsen, N., 1997b. Geomagnetic Tides and Related phenomena, in, Lecture
notes in Earth Sciences, 66, ed. Wilhelm, H., Springer-Verlag, Berlin.
Olsen, N., 1998. The electrical conductivity of the mantle beneath Europe
derived from C-Responses from 3 h to 720 h, Geophys. J., 133, 298
308.
Olsen, N., 1999. Induction studies with satellite data, Surv. Geophys., 20,
309340.
Olsen, N. et al. 2000. rsted initial eld model, Geophys. Res. Lett., 27,
36073610.
Olson, W.P., 1984. Introduction to the topology of magnetospheric current
systems, in Magnetospheric Currents, pp. 4962, ed. Potemra, T.A., Am.
geophys. Un. Monogr., 28.
Onwumechilli, C.A. & Ozoemena, P.C., 1989. Contours of equatorial elec-
trojet current density, J. Atmos. Terr. Phys., 51, 163168.
Parkinson, W.D. & Hutton, V.R.S., 1989. The electrical conductivity of the
Earth, in Geomagnetism, Vol. 3, pp. 261322, ed. Jacobs, J.A., Academic
Press, London.
C
2002 RAS, GJI, 151, 3268
68 T. J. Sabaka, N. Olsen and R. A. Langel
Peredo, M., Stern, D.P. & Tsyganenko, N.A., 1993. Are existing magne-
tospheric models excessively stretched?, J. geophys. Res., 98, 15 343
15 354.
Potemra, T.A., 1982. Birkeland currents: present understanding and some re-
maining questions, in High-latitude Space Plasma Physics, eds Hultqvist,
B. & Hagfors, T., Nobel Symposium, 54, Plenum, New York.
Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., 1992. Nu-
merical recipes in FORTRAN: the Art of Scientic Computing, 2nd edn,
pp. 1963, Cambridge University Press, Cambridge.
Purucker, M.E., Sabaka, T.J., Langel, R.A. & Olsen, N., 1997. The missing
dimension in Magsat and POGO anomaly studies, Geophys. Res. Lett.,
24, 29092912.
Rangarajan, G.K., 1989. Indices of geomagnetic activity, in Geomagnetism,
Vol. 3, pp. 323384, ed. Jacobs, J.A., Academic Press, London.
Ravat, D.N., Langel, R.A., Purucker, M.E., Arkani-Hamed, J. & Alsdorf,
D.E., 1995. Global vector and scalar Magsat magnetic anomaly maps,
J. geophys. Res., 100, 20 11120 135.
Richmond, A.D., 1995. Ionospheric electrodynamics using magnetic apex
coordinates, J. Geomag. Geoelectr., 47, 191212.
Sabaka, T.J. & Baldwin, R.T., 1993. Modeling the Sq magnetic eld from
POGO and Magsat satellite and cotemporaneous hourly observatory
data: phase I, Contract Report HSTX/G&G9302, Hughes STX Corp.
for NASA/GSFC Contract NAS531 760.
Sabaka, T.J., Langel, R.A., Baldwin, R.T. & Conrad, J.A., 1997. The geo-
magnetic eld 19001995, including the large-scale eld from magneto-
spheric sources, and the NASA candidate models for the 1995 revision of
the IGRF, J. Geomag. Geoelectr., 49, 157206.
Sabaka, T.J., Olsen, N. & Langel, R.A., 2000. A comprehensive model
of the near-Earth magnetic eld: phase 3, NASA/TM-2000-209894, 1
75.
Schmitz, D.R., Meyer, J. & Cain, J.C., 1989. Modeling the Earths geomag-
netic eld to very high degree and order, Geophys. J., 97, 421430.
Schmucker, U., 1985. Magnetic and electric elds due to electromagnetic
induction by external sources, in Landolt-B ornstein, New-Series, 5/2b,
Springer-Verlag, Berlin.
Schmucker, U., 1999. Aspherical harmonic analysis of solar daily variations
in the years 196465-I. Methods, Geophys. J. Int., 136, 439454.
Schumaker, L.L., 1981. Spline Functions: Basic Theory, Wiley, New York.
Seber, G.A.F. & Wild, C.J., 1989. Nonlinear Regression, pp. 1768, Wiley,
New York.
Takeda, M. & Maeda, H., 1983. F-region dynamo in the evening
interpretation of equatorial LD anomaly found by Magsat, J. Atmos. Terr.
Phys., 45, 401408.
Tarantola, A. & Valette, B., 1982. Generalized non-linear inverse problems
solved using the least squares criterion, Rev. Geophys., 20, 219232.
Toutenburg, H., 1982. Prior Information in Linear Models, pp. 1215, Wiley,
New York.
Tsyganenko, N.A., 1987. Global quantitative models of the geomagnetic
eld in the cislunar magnetosphere for different disturbance levels, Planet.
Space Sci., 35, 13471358.
Tsyganenko, N.A., 1989. A magnetospheric magnetic eld model with a
warped tail current sheet, Planet. Space Sci., 37, 520.
Volland, H., 1984. Atmospheric Electrodynamics, pp. 1205, Springer-
Verlag, Berlin.
Walker, A.D. &Backus, G.E., 1997. Asix-parameter statistical model of the
Earths magnetic eld, Geophys. J. Int., 130, 693700.
Wessel, P. & Smith, W.H.F., 1991. Free software helps map and display data,
EOS, Trans. Am. geophys. Un., 72, 441.
Winch, D.E., 1981. Spherical harmonic analysis of geomagnetic tides, 1964
1965, Phil. Trans. Roy. Soc. Lond., A., 303, 1104.
C
2002 RAS, GJI, 151, 3268

You might also like