Data-based modelling of the Earth ’ s dynamic magnetosphere : a review

This paper reviews the main advances in the area of data-based modelling of the Earth’s distant magnetic field achieved during the last two decades. The essence and the principal goal of the approach is to extract maximum information from available data, using physically realistic and flexible mathematical structures, parameterized by the most relevant and routinely accessible observables. Accordingly, the paper concentrates on three aspects of the modelling: (i) mathematical methods to develop a computational “skeleton” of a model, (ii) spacecraft databases, and (iii) parameterization of the magnetospheric models by the solar wind drivers and/or ground-based indices. The review is followed by a discussion of the main issues concerning further progress in the area, in particular, methods to assess the models’ performance and the accuracy of the field line mapping. The material presented in the paper is organized along the lines of the author Julius-Bartels’ Medal Lecture during the General Assembly 2013 of the European Geosciences Union.


Introduction
The geomagnetic field is the principal agent connecting our planet's ionosphere with the highly variable interplanetary medium, incessantly disturbed by dynamical processes at the Sun.The Earth's magnetosphere serves as a giant storage reservoir of energy pumped in from the solar wind and intermittently spilled into the upper atmosphere during space storms.As humankind becomes more and more dependent on space technologies, it becomes increasingly important to be able to accurately map the distant geomagnetic field and predict its dynamics using data from upstream solar wind monitors.Two approaches to the problem have been successfully pursued over recent decades.The first is to treat the solar wind as a flow of magnetized conducting fluid and to numerically solve first-principle equations, governing its interaction with the terrestrial magnetic dipole.Based on pure theory, that approach addresses the question: "What would the magnetosphere look like and how would it behave if the underlying approximations and techniques were universally accurate?"This review focuses on the other, completely different approach, based on direct observations.Its essence is to develop an empirical description of the global geomagnetic field and its response to solar wind driving by fitting model parameters to large multi-year sets of spacecraft data.Models of that kind seek to answer the question: "What can in situ measurements tell us about the global magnetospheric configuration and its storm-time dynamics, provided our approximations are realistic, flexible, and the data coverage is sufficiently dense and broad?"Five decades of spaceflight have produced enormous amounts of archived data and a whole suite of empirical models have already been developed on that basis (e.g., McCollough et al., 2008, and references therein).Recent and ongoing multi-spacecraft missions keep pouring in new data and further expand the huge and yet largely untapped resource of valuable information.The main goal of such data-based modelling is to extract the largest possible knowledge from the accumulated data, thus synergistically maximizing the output of present and past space experiments.Most of the existing models of this kind are implemented as self-contained computer codes, available to N. A. Tsyganenko: Data-based models of the magnetosphere the magnetospheric community as relatively simple handson tools for researchers, which have proved to be a useful by-product of our efforts.
This paper presents a condensed overview of methods and results of empirical magnetosphere modelling.It can be viewed as an update on the previously published review articles on the subject (Tsyganenko, 1990;Stern, 1994).The last two decades have seen significant advances in the field, such that the above cited reviews have become largely obsolete.The purpose of the present paper is to highlight novel techniques, summarize recent progress, and outline basic directions for future research.
The article is organized as follows.Section 2 briefly discusses the main differences between the modelling of the internal and external components of the total magnetospheric magnetic field.Section 3 is devoted to the mathematical structure of existing models and outlines basic methods to represent contributions to the observed field from the principal magnetospheric field sources.Section 4 addresses methods to parameterize the empirical models, that is, to relate the magnitudes and geometrical characteristics of the field sources to routinely monitored external input variables.Section 5 is a short overview of spacecraft data contained in available archives, basic preparation procedures, and requirements to be met for the data to be included in the modelling sets.Section 6 is focused on issues of the models' performance and accuracy.Section 7 examines the problem of consistency between the empirical model B field and distributions of the magnetospheric plasma pressure.Finally Sect.8 discusses outstanding problems and challenges to be addressed in future data-based modelling studies.

Internal and external parts of the magnetospheric magnetic field
The total magnetospheric magnetic field vector B can be represented as the sum of the internal part B I (also called the "main geomagnetic field"), and the external part B E , associated with electric currents flowing inside and outside Earth, respectively.Note that, following the insightful analysis by Vasyliunas (2005), we intentionally used above the term "associated" instead of "produced" (or "generated"), keeping in mind that, strictly speaking, in astrophysical objects like the magnetosphere, electric currents should be viewed as a result of the interaction between the bulk of solar wind plasma and the magnetospheric magnetic field, rather than its source.Nonetheless, for the sake of brevity, in the following we will retain the short term "sources" for the magnetospheric currents, tacitly keeping in mind its conditional and relative meaning.for the scalar potential U (r, θ, φ) of the main geomagnetic field has remained virtually unchanged since then, except for its length: modern IGRF models include the terms up to 10th order in n, owing to the dramatically increased flow of data from a large number of ground-based observatories, complemented by a huge volume of marine, airborne, and satellite data.For comprehensive information on the main field modelling, we refer the reader to topical reviews (e.g., Langel, 1987) and monographs (Chapman and Bartels, 1962).
The internal part B I = −∇U largely dominates on the ground and at low altitudes, but rapidly decreases with geocentric distance as ∼ r −3 and becomes comparable to the external field at r ∼ 10 R E (on the order of the magnetopause standoff position).Beyond that distance (in the magnetotail), the external part comes into foreground, while the internal part asymptotically falls off to zero.
Due to its curl-free nature above the Earth's surface, the internal field B I is uniquely defined in the entire space by the expansion coefficients in Eq. ( 2), which can be accurately computed using only ground-based and low-altitude data.By contrast, the external field B E is associated with volume currents widely distributed over the magnetosphere, and hence cannot be described by a scalar potential.This means that, unlike the internal part, the global external field cannot be derived from spatially localized observations, and this is why empirical magnetosphere models critically depend on extensive sets of spacecraft data, covering the modelling region with a sufficient density.
Another important complication is that, unlike the internal field which is almost static in the Earth's frame of reference (barring slow secular variations), the external field is highly variable over a wide range of timescales -from seconds, to hours, days, and up to the 11 yr solar cycle period.The variations are due to several factors, such as the Earth's rotation and its orbital motion around Sun, resulting in diurnal and yearly oscillations of the geodipole orientation with respect to the Sun-Earth line, incessant changes in the state of the incoming solar wind flow, and irregular internal instabilities in the magnetosphere.
Last but not least is the fundamental difference in the nature of data sets used in the modelling of B I and B E .In the former case, simultaneous data from a host of ground and low-altitude locations are available in almost real time.This makes it possible not only to create and periodically update accurate models of the main geomagnetic field, but also to dynamically reconstruct ground variations of the field of low-altitude sources, such as the ionospheric and fieldaligned currents (henceforth, FACs, for short).By contrast, in the case of the distant B E the situation is quite the opposite: the highly variable field occupying a huge domain can be measured, in the best case, at only a few locations at a time.Regretfully, the project to simultaneously monitor the magnetosphere by a widely distributed swarm of 50-100 space probes (Angelopoulos et al., 1998) still remains in the realm of dreams.The principal goal of empirical modelling is to partially overcome this difficulty, by taking advantage of the abundance of archived space magnetometer and plasma instrument data from many past and ongoing missions, covering a wide variety of diverse magnetospheric events.

Mathematical framework of data-based models
If one likens empirical models to a building structure, then it can be said to rest on three pillars.The first pillar is the mathematical framework, i.e., a set of equations representing contributions to the total field of individual magnetospheric current systems.The second pillar is the spacecraft and groundbased data, used to determine optimal values of model parameters.The third pillar is the parameterization methods and equations, relating the magnitudes and geometrical characteristics of individual field sources, as well as their temporal dynamics, to routinely available parameters of the incoming solar wind and/or ground geomagnetic activity indices.
This section outlines basic principles and methods to mathematically represent contributions to the external field from individual magnetospheric current systems.Most of the following material corresponds to the advanced approach that has been developed in the past decade (e.g., Tsyganenko, 2002a, b;Tsyganenko andSitnov, 2005, 2007, andreferences therein).It should be noted from the outset that, from the viewpoint of physics, magnetospheric currents actually form a single entity.Dividing them into separate components is largely a matter of convenience, justified by the fact that different parts of the whole current system have different geometry, differently respond to external driving, and have largely different relaxation timescales.It has been commonly accepted to represent the net external field B E as a sum of contributions from the ring current, B RC , tail current sheet, B TC , large-scale field-aligned current systems, B FAC (including both Region 1 and 2), and the magnetopause currents, B MP , so that the total field (3) Note that in all recent models (T96 and later) the above expansion also included the so-called "interconnection" field B INT , proportional to the transverse component of the IMF.Adding that term was motivated by the well-known fact that the IMF partially penetrates into the magnetosphere, most conspicuously manifested in the correlation of the B y field components (Fairfield, 1979;Cowley, 1981;Cowley and Hughes, 1983;Sergeev, 1987).This question will be further discussed in more detail in Sect.8.1.
The magnetopause field B MP is not an independent term: it is added to all other parts of the total B vector to ensure full confinement (or "shielding") of the magnetospheric magnetic field inside the common model boundary S, so that where n is unit normal vector to the magnetopause.Starting from the T96 model (Tsyganenko, 1995(Tsyganenko, , 1996)), and in all later data-based models, the magnetopause S has been represented by an independently pre-defined empirical surface, fitted to data of boundary crossings by satellites, which makes the boundary condition Eq. ( 4) linear with respect to B. This prompts us to split the term B MP into a sum of partial fields, each of which serves as a shielding field for the corresponding term (of the first four) in the right-hand side of Eq. ( 3), so that the total field reads where each of the four paired terms is independently shielded within the boundary.As detailed in the following sections, a natural way to increase the model's flexibility is to further expand the partial fields B RC , B TC , and B FAC , representing them as linear combinations of independent normalized vector fields b TC , and b FAC , paired with their respective shielding fields h

TC , and h (k)
FAC .As a result, in the most general case the field of each (ith) source assumes the generic form of an expansion where each kth term in the sum includes a linear coefficient a (k) i and a set of nonlinear parameters {α (k) i }, quantifying the magnitude and geometrical properties of the partial field source, as well as its response to the model's input quantities, including the geodipole tilt angle , the solar wind speed and dynamic pressure P dyn , the interplanetary magnetic field (IMF), and related external driving variables.Each term in Eq. ( 6) satisfies the shielding condition at the magnetopause S b which is the principal advantage of the approach, since it makes it possible to independently vary the parameters of individual magnetospheric field sources and, at the same time, keep the total field fully shielded inside S for any values of the coefficients {a (k) i } and (within a certain finite range) of the variable nonlinear parameters {α The first pair of terms in Eq. ( 5), corresponding to the shielded Earth's main field, is treated separately.The internal field B I is known in advance with great accuracy from IGRF expansions, and, once a model magnetopause shape and size is known, the corresponding shielding field B MP,I can be uniquely obtained in a straightforward way.Since the magnetopause is located relatively far from Earth, all higherorder harmonics of the main field are small there, so that B I can be accurately approximated by a purely dipolar field and, hence, the only quantities that control B MP,I are the dipole tilt angle and the solar wind parameters that define the size and shape of the model boundary.In Sect.3.4 we will address the derivation of the shielding fields in greater detail.

Equatorial magnetospheric currents and their magnetic field
From a global viewpoint, the observed magnetospheric B field structure is shaped by two plasma domains: (i) the magnetosheath and the polar cusps (which themselves can be viewed as extensions of the magnetosheath inside the dayside magnetosphere), and (ii) the nightside equatorial region, from the outer boundary of the inner magnetosphere to the distant tail plasma sheet.In the empirical approach to magnetic field modelling we disregard the issue of consistency between the magnetic field and plasma pressure (that subject is addressed in more detail in Sect.7) and represent the model field by a formal superposition of analytically simple modules.
Physically, the inner ring current and the more distant tail current sheet form a single equatorial current.In a topological sense, the difference between the two is that the ring current flow lines encircle Earth and are fully closed inside the magnetosphere, whereas the tail currents flow in the azimuthal direction within a limited sector of longitudes and then close via the magnetopause, forming "theta"shaped current loops.Nevertheless, when constructing a fully shielded magnetic field model, both the ring and tail currents can be regarded as laterally unbounded equatorial sources, extending arbitrarily far beyond the magnetopause.This is illustrated in Fig. 1, where the top pair of panels show spatially unrestrained electric current flow lines (red traces in the 3-D view on the left) and corresponding lines of the unshielded magnetic field B TC in the noonmidnight meridian plane (blue traces on the right), extending beyond the model magnetopause (grey-shaded surface and purple line).Adding the field B MP,TC results in full confinement of the shielded field within the magnetopause, so that the total normal component B TC + B MP,TC • n S = 0 everywhere on the boundary.Now the magnetic field (hence the electric currents) outside the magnetosphere can be nullified without violating Maxwell's equations; the resulting jump in the previously continuous tangential field component will correspond to a surface current, exactly equal to that needed to redirect the equatorial current and close it over the boundary, as illustrated in the bottom panels of Fig. 1.The above described "gedanken experiment" was first realized by Stern (1987, Appendix A) and further substantiated by Sotirelis et al. (1994).

Modelling the ring current field
The ring current is a principal source of the external field in the inner magnetosphere, in particular during storms when it dramatically grows in magnitude and becomes strongly asymmetric due to the formation of a duskside partial ring current.In early empirical models (Tsyganenko and Usmanov, 1982;Tsyganenko, 1987Tsyganenko, , 1989; henceforth, TU82, T87, and T89) the ring current field was represented by a very compact two-parameter axisymmetric module, based on a simple modification of the dipolar vector potential, expressed in cylindrical coordinates {ρ, φ, Z} as A = Ae φ , with A = 4B 0 ρ 3 0 ρ(Z 2 + ρ 2 + 4ρ 2 0 ) −3/2 .The model was parameterized by the scale radius ρ 0 and the scale intensity B 0 , equal to the model field magnitude at the origin.In the later T96 model, both the ring and tail current fields were represented by more sophisticated potentials (see Sect. 3.1.2below), arranged in combinations of several terms in order to confine the currents within a limited range of radial distance and the Z coordinate.
The above-referenced solutions can be used as building blocks in constructing more realistic fields, taking into account, for example, the eastward current due to the positive radial gradient of the particle pressure in the innermost region at r ≤ 2 − 3 R E .Unfortunately, all these models are axially symmetric, while, as already said, the actual ring current can develop a strong asymmetry during storms.The azimuthal asymmetry of the particle pressure results in the divergence of the equatorial current and formation of fieldaligned, or Birkeland, currents.As a result, the problem becomes three-dimensional, and to devise a realistic solution we need to turn to theory.
Since the ring current flows relatively close to Earth, where the total magnetic field is not drastically different from its main (dipolar) component, one can calculate the drift and magnetization electric current densities j d and j m as where the perpendicular and parallel particle pressures P ⊥ (r e , φ) and P (r e , φ) are a priori defined as functions of the equatorial radial distance r e and the longitude φ.Note that, strictly speaking, the pressures P ⊥ (r e , φ) and P (r e , φ) and the "background" magnetic field B o should be mutually consistent, in other words, must form a force-balanced configuration.Nevertheless, in the low-beta approximation, one still can use Eqs.( 8) and ( 9) to roughly calculate the currents in an apriori prescribed magnetic field.Thus obtained currents are then used to evaluate the associated disturbance magnetic field.That problem was addressed in many works, starting from the pioneering study by Akasofu and Chapman (1961), and followed by successful attempts to iteratively derive higher-order solutions, taking into account the perturbation field of the ring current itself (e.g., Sckopke, 1972, and references therein).All those studies used a purely dipolar background field as a starting approximation for the background field B o , and employed the above gyrotropic equations ( 8)-( 9) for the electric currents.A notable exception in this sense was a work by Lackner (1970), based on a more general Vlasov formalism.
The first problem with the above models is that they were limited to axially symmetric plasma configurations with ∂P ⊥ /∂φ = ∂P /∂φ = 0 and, for that reason, they did not include FACs.The FACs can be evaluated (e.g., Birmingham, 1992a, b) by integrating the divergence of the drift current along a field line connecting the point s, where the current j is to be calculated, with a magnetically conjugate location e in the equatorial plane However -and this is the second problem -for the purposes of data-based field modelling, it is not enough to simply numerically evaluate the magnetic field of the ring current.This is only the first step, while the greatest challenge and the final goal is to obtain a reasonably compact and flexible global analytical description of the disturbance field, which can be fitted to satellite data.Both the above issues were first addressed in (Tsyganenko, 2000), where azimuthally asymmetric particle pressure distributions were used to calculate the first-order drift, magnetization, and field-aligned currents.The essence of the approach was to separately represent the symmetric and partial components of the ring current, by specifying for each part its own distribution of the equatorial plasma pressure.The symmetric ring current (SRC) was treated as a basic permanent feature of the inner magnetosphere, and the corresponding radial distribution of the plasma pressure was assumed in the form of smooth analytical approximations for P (SRC) ⊥ (r e ) and for the anisotropy parameter γ (r e ) = P (SRC) ⊥ /P (SRC) .Both profiles were fitted by least squares to quiet-time experimental curves by Lui and Hamilton (1992), in which the pressure peaks at r e ∼ 2.8 R E .
Storm-time variations were supposed to be reproduced by varying the magnitude and scaling the size of the SRC.Unlike the SRC, the partial ring current (PRC) develops to its full extent only during active periods, owing to enhanced plasma convection from the tail.For that reason, the PRCrelated pressure P (PRC) was assumed to be isotropic and peaked at larger distances, around r e ∼ 6-7 R E .Its variation with longitude φ was represented by a sum of lowest-order Fourier terms, so that where the radial variation is factored out in P (PRC) 0 (r e ), the parameter ε controls the degree of azimuthal asymmetry, and the phase angle φ 0 defines the longitude of the PRC peak. Figure 2 illustrates the configuration of electric current flow lines, obtained from Eqs. ( 8)-( 11) as a superposition of the axisymmetric and "quadrupole" PRC components, corresponding to the first and second bracketed terms in the pressure Eq. ( 11).
The current densities were calculated using a purely dipolar background magnetic field B o , which eliminated the need to numerically trace the field lines in the calculation of the electric currents from Eqs. ( 8)-(11).In addition, the axial symmetry of the dipolar B o , combined with the purely harmonic azimuthal variation of the pressure in Eq. (11) made it possible to reduce the problem to 2-D.These two factors allowed us to represent the SRC and PRC fields using computationally fast analytical approximations, included later on in the T02 (Tsyganenko, 2002a, b) and TS05 (Tsyganenko and Sitnov, 2005) empirical models.Their relative simplicity, however, came not without a price: using a purely dipolar background field resulted in inaccurate mapping between the equatorial PRC and Region 2 (R2) FACs at low altitudes.An advanced PRC model based on a realistic asymmetric background field (Tsyganenko, 2013, referred to henceforth as T13) yields more accurate results, but demands much more computing resources.

Modelling the magnetic field of the tail current
There exists a wide variety of analytically simple magnetic fields associated with planar current sheets and disks.One can start, for example, from the simplest source in the form of a straight linear current, flowing in the equatorial plane parallel to the Y GSM axis at X = X 0 , which spreads out in (re), the parameter ε controls the degree of azimuthal asymmetry, and the phase angle φ0 defines the longitude of the PRC peak. Figure 2 illustrates the configuration of electric current flow lines, obtained from Eqs.( 8)-( 11) as a superposition of the axisymmetric and 'quadrupole' PRC components, corresponding to the first and second bracketed terms in the pressure equation ( 11).The current densities were calculated using a purely space over a scale half-thickness D. Its field can be represented by the elementary vector potential dA = dA y e y , where dA y Integrating it over X 0 with different weight functions I (X 0 ) provides a family of simple analytical fields, corresponding to spread-out current sheets with a finite half-thickness D, with various radial profiles of the electric current density I (X).In the TU82 model, a linear variation of I (X 0 ) was assumed between the inner and outer edges of a planar current sheet, which yielded a simple magnetotail field module.A more sophisticated hyperbolic form of I (X 0 ) was adopted in the T87 model, which made it possible to extend its validity range further out into the distant tail.
Several other simple functions I (X 0 ) can be found, which yield the corresponding magnetic field components in a closed form.One example is a bell-shaped current density profile, centered at which results in a compact vector potential with only an A ycomponent in the form Dividing Eq. ( 13) by X and differentiating the result with respect to that parameter yields another solution, which differs from the original one by much steeper slopes of the bellshaped profile.Such a current "slab" module was used in the T13 model to improve its flexibility in the dayside sector.Note that the parameter D in Eq. ( 13) can be assumed to be a function of coordinates, making it possible to model spatial variations of the current sheet thickness.Another family of remarkably simple analytic solutions for the magnetic field, widely used in empirical modelling, is associated with axially symmetric disk-like equatorial distributions of the electric current (Tsyganenko, 1989(Tsyganenko, , 1990)).It is derived by equating to zero the electric current density outside an infinitely thin current sheet, expressed in cylindrical coordinates {ρ, φ, Z} via the azimuthal component of the vector potential A = A(ρ, Z)e φ .A general solution of the corresponding 2nd-order equation ∇ × (∇ × A) = 0 reads from which the weight function C(K) is derived by applying Bessel's transform to the B z -component of the equatorial field, corresponding to the potential Eq. ( 14).Specifying B z (ρ) as a simple bell-shaped profile of the magnetic field depression centred at the origin, B z (ρ) ∼ (ρ 2 +a 2 ) −1/2 , leads to a compact solution for the potential where the parameter a defines a characteristic scale length of the current density radial profile, and Due to the presence of |Z|, the above potential exhibits a kink at the plane Z = 0, corresponding to infinitely thin current sheet.Replacing |Z| by ζ = √ Z 2 + D 2 spreads the thin sheet over a finite bell-shaped profile with a scale halfthickness D, which can be further made a function of coordinates, allowing one to model magnetic fields of current disks with a variable thickness.Successive differentiation of Eq. ( 15) with respect to a yields a sequence of independent vector potentials with progressively faster rates of asymptotic decrease of the current with growing radial distance.Final equations for the first three potentials A (1) , A (2) , and A (3) are To save page space, we omit the corresponding equations for the field components, which can be easily derived by calculating ∇ × A. The set of solutions described above can be either directly used to generate independent modules b (k) i in Eq. ( 6), or can be first arranged into linear combinations with the coefficients and scale lengths defined in such a way that they form a set of ad hoc modules with desirable radial profiles of the electric current.The latter approach was adopted in the T96, T02, and TS05 models, though using somewhat different basic potentials.
The rapidly growing volume of archived space magnetometer data suggests the need to look for ways to enhance the models' capability to ingest new information and reproduce the structure of the magnetosphere in more detail.In the modelling of the main geomagnetic field, this can be done simply by adding more higher-order harmonics into the scalar potential expansion (2).An interesting and important question is whether a similar approach could be developed and implemented in the external field modelling.The first step in that direction was made by Tsyganenko and Sitnov (2007), who devised the TS07D model, based on extensible high-resolution expansions for the field of the equatorial current sheet.The original idea was to start from vector potentials in the integral form Eq. ( 14), but instead of transforming them to particular closed-form solutions like Eq. ( 16), replace the integrals by formal expansions over a discrete equidistant set of wavenumbers K, in which higher values of K would correspond to smaller-scale details in the current sheet structure.A particular problem, however, was that the above single-component vector potential A = A(ρ, Z)e φ represents only the axisymmetric part of the field, corresponding to purely azimuthal equatorial currents, whereas in the general case the model must also include azimuthally asymmetric terms, which is especially important to accurately describe a pronounced duskside depression of the storm-time field in the inner magnetosphere.
Mathematically, introducing azimuthally asymmetric terms is not quite straightforward, since in that case both the electric current and the vector potential can no longer be described by a single-component vector.The problem was circumvented by starting from the outset with scalar potentials and then converting them into vector potentials.More specifically, we first solve Laplace's equation for magnetic scalar potentials northward and southward from an infinitely thin equatorial current sheet, then transform the obtained scalar potentials γ + and γ − into a single vector potential, and finally, modify it by spreading the originally infinitely thin current sheet over a finite thickness across the equatorial plane.Details of the derivation can be found in the original paper; here we reproduce only the final form of the expansion terms.The axially asymmetric terms include both factors sin(mφ) and cos(mφ) with m = 1, 2, . . ., responsible for the noon-midnight and dawn-dusk asymmetries, respectively, as follows: The axisymmetric term is represented separately as Replacing the integration in Eq. ( 14) by a summation over a discrete spectrum of the wavenumbers k n yields the model expansion for the potential The method outlined above makes it possible to reveal some interesting details of the storm-time dynamics of the magnetospheric currents.Those results will be discussed in Sect.4.1.2below.Here we only note that, unlike all the earlier models with "custom-tailored" ring current and tail field modules, individual terms in the expansion Eq. ( 19) should be viewed simply as formal Fourier terms, and hence cannot be associated with any specific mode of the external driving and internal decay, as was the case in the TS05 model.A more detailed discussion of this issue is deferred to Sect.4.1.2.

Modelling the magnetic field associated with large-scale Birkeland currents
In this section we address mainly the effects of the Region 1 (R1) field-aligned currents, since the R2 currents should be viewed as an intrinsic part of the PRC, already discussed in Sect.3.1.1.The R1 currents form the outermost internal current system in the magnetosphere, topologically closest to the magnetopause, and serve as the shortest link between the solar wind generator in the magnetosheath and the high-latitude ionosphere.The greatest problem in the empirical modelling of the R1 FACs is that there exists no satisfactory quantitative theory of those currents which could elucidate their geometry in the distant magnetosphere.While the R2 currents and their magnetic field can be modeled using static force balance equations and the observed pressure distributions in the closed field line region, the only way to represent the effects of the R1 currents is to empirically specify a flexible model and derive its parameters from available data.
A convenient way to define a current system geometry is to represent the corresponding volume density via the Euler potentials which automatically guarantees its continuity, ∇ • j = 0.As concerns the R1 FACs, the only piece of experimental evidence we can rely upon is that at low altitudes they flow into and out of the ionosphere along quasi-dipolar field lines, and their intersection with the ionosphere is an eccentric band, which roughly matches the auroral oval.Details of the ionospheric closure of the FACs are of little interest in our case, because of the negligible effect of the ionospheric currents beyond R ≥ 1.5 − 2 R E (Tsyganenko, 2002a).Based on the above, the Euler potentials in Eq. ( 20) should be defined in such a way that the electric current flow lines nearly follow the dipolar magnetic field lines at low altitudes, but deviate from them at larger distances.Since we have no a priori knowledge on the FAC geometry in the distant magnetosphere and intend to extract that information from data, the potentials ξ and χ must be sufficiently flexible.Suitable functions satisfying the above requirements were introduced by Tsyganenko and Stern (1996) in solar-magnetic spherical coordinates {r, θ, φ} as where and the second term i (φ) in Eq. ( 21) is the colatitude of the R1 oval at ionospheric altitude, defined as a periodic function of longitude φ.A fundamental property of the surface ξ(r, θ, φ) = 0, which can be easily verified from Eqs. ( 21)-( 22), is that at low altitudes its shape is close to the surface formed by dipolar field lines, crossing the Earth's surface along the oval = i (φ), while at large radial distances it asymptotically approaches the equatorial plane.The parameter α defines the location of the transition region between the dipole-like and tail-like shape of meridional cross sections of the surface.Larger values of α correspond to a larger curvature of the surface and vice versa, as shown in Fig. 3, which displays meridional sections of surfaces of constant (r, θ ) for two values of α.
The second potential χ(r, φ) defines the shape of the electric current flow lines on the surface ξ(r, θ, φ) = 0 as well as the azimuthal distribution of the FAC density.Radially independent potentials with ∂χ/∂r = 0 correspond to purely poloidal currents, flowing in the meridional planes φ = const.Introducing a radial variation in χ adds an azimuthal component to j , which can be used to make the nightside FACs either exit the magnetosphere via its flanks or close across the midnight meridian.
The method outlined above was used to construct a flexible R1 FAC module in the T96 model (the first one to explicitly include the FAC contribution), as well as in more recent T02, TS05, and TS07 models.In the T96 model, the second potential was assumed in the form where the radial G(r) and azimuthal f (φ) factors were a priori defined to make the distant FACs enter and exit the magnetosphere via its the dawnward and duskward flanks, and to place the duskside and dawnside peaks of the lowaltitude FACs closer to the noon meridian.Upon having defined the spatial distribution of FACs, their magnetic field was computed by Biot-Savart integration in the entire modeled region of the magnetosphere, then individually approximated by suitable potential fields in the high-latitude and low-latitude current-free domains, and finally, interpolated across the transition regions (FAC sheets) separating these domains.
In the T02, TS05, and TS07 models (see Tsyganenko, 2002a, for details), the procedure was somewhat different.Instead of specifying from the outset the azimuthally asymmetric FAC sheet, we started from the axially symmetric surface (r, θ ) = i0 with the constant ionospheric colatitude i0 of the R1 zone.The second potential χ was assumed to be a function only of the longitude φ in the simplest form f (φ) = sin mφ, with the goal to represent the local time distribution of the FACs by the first few Fourier harmonics.As in the case of the PRC (see Sect. 3.1.1above) the assumed axial symmetry of the surface and the sinusoidal variation of the FAC density greatly simplified the problem by making it possible to isolate the φ-dependence in the magnetic versa, as shown in Fig. 3, which displays meridional sections of surfaces of constant Θ(r, θ) for two values of α.The second potential χ(r, φ) defines the shape of the electric current flow lines on the The method outlined above was used to construct a flexible R1 FAC module in the T96 model (the first one to explicitly include the FAC contribution), as well as in more recent T02, TS05, and TS07 models.In the T96 model, the second potential was assumed in the form where the radial G(r) and azimuthal f (φ) factors were a priori defined to make the distant FACs enter and exit the magnetosphere via its the dawnward and duskward flanks, and to place the duskside and dawnside peaks of the low-altitude FACs closer to the noon meridian.Upon having defined the spatial distribution of FACs, their magnetic field was computed by Biot-Savart integration in the 420 entire modeled region of the magnetosphere, then individually approximated by suitable potential fields in the high-latitude and low-latitude current-free domains, and finally, interpolated across the transition regions (FAC sheets) separating these domains.
14 field components into separate factors sin mφ and cos mφ (Tsyganenko, 1993, Appendix B) which allowed us to reduce the problem to 2-D and define the entire 3-D field by calculating its components by Biot-Savart integration in only a single meridian plane.The next step was to analytically approximate the obtained field, which was done by starting from a family of so-called "conical" harmonics, representing the field of a conical current sheet (Tsyganenko, 1991), and modifying that field by a suitable 2-D deformation in spherical coordinates (the essence of the deformation method is described below in Sect.3.3).The final step was a rotational deformation around the y axis, introduced to replicate the observed day-night asymmetry of the global system of R1 FACs. Figure 4 illustrates the geometry of the electric current flow lines, obtained from the shielded R1 FAC model field by numerically calculating its curl.As verified by calculations of the R1 FAC dynamics in real events (see an example in Sect.8.3, Fig. 19), the model outlined above yielded quite reasonable results, and was implemented in the form of relatively fast numerical codes.However, the approach is not free from drawbacks.For one thing, the model is insufficiently flexible; in particular, varying the parameter α in Eq. ( 22) could help to derive from data the optimal shape of the R1 FAC surface.Unfortunately, in the framework of the T02/TS05 approach it would require one to iteratively recalculate the entire set of deformation parameters, a computationally unfeasible task.
Another deficiency of the deformation method is that it works fairly well only for the lowest 1st and 2nd Fourier harmonics, but rapidly deteriorates for higher-order terms, which does not allow one to model the magnetic effects of azimuthally localized Birkeland currents.Also, it is a priori assumed in the model that the currents have no azimuthal component, which prevents one to explore the FAC closure in the distant magnetosphere.
An alternative way to build a model with much greater flexibility is to evaluate the magnetic field B FAC due to FACs using Biot-Savart integration.In terms of a vector potential, the problem reduces to calculating In the T02, TS05, and TS07 models (see (Tsyganenko, 2002a)  in only a single meridian plane.The next step was to analytically approximate the obtained field, which was done by starting from a family of so-called 'conical' harmonics, representing the field of a conical current sheet (Tsyganenko, 1991), and modifying that field by a suitable 2D deformation in spherical coordinates (the essence of the deformation method is described below in Section 3.3).The final step was a rotational deformation around the Y-axis, introduced to replicate the observed daynight asymmetry of the global system of R1 FACs.
over a set of electric current flow lines, aligned with the background magnetic field B 0 .Using the vector potential makes it possible to conserve ∇ • B = 0 and regularize the integrand in Eq. ( 24) by introducing a finite transverse scale length D = D(s), so that the denominator takes the form (Tsyganenko, 1997(Tsyganenko, , 2000)).To keep the currents field-aligned, it suffices to set D(s) ∼ B −1/2 0 , so that the magnetic flux inside the electric current flow tube remains constant.To speed up the computation of the integral in Eq. ( 24), the multitude of volume elements constituting the smooth FAC flow tube can be replaced by much smaller sets of straight segments with linearly varying half-thickness D(s), with the vector potential components expressed in a closed analytical form.Making the segment lengths proportional to the local curvature radius of the current flow tube dramatically decreases the number of summation terms, as sketched in Fig. 5.We used that method in constructing a numerical model of the substorm current wedge (Sergeev et al., 2011).A similar approach was developed independently by Ontiveros et al. (2006) in their model of the R2 FACs and the magnetopause currents.As an illustration, Fig. 6 presents a sample test distribution of the model R1 and R2 FAC density at the ionospheric level, obtained using the above summation procedure.A great advantage of the Biot-Savart summation is that, unlike in the magnetic field deformation method, it allows one to modify the current system geometry without adding unwanted artificial currents.On the negative side, in most cases the procedure is computationally rather intensive, even in its finite-segment version.

Geodipole tilt effects
The Earth's dipole axis is inclined by ≈ 10 • (as of the 2010 epoch) to its rotation axis, which, in its turn, is inclined by 23.4 • to the normal to the ecliptic plane.This results in diurnal and yearly variations of the angle between the geodipole axis and the terminator plane within the range −33.4 • ≤ ≤ 33.4 • .The dipole tilt variations affect help to derive from data the optimal shape of the R1 FAC surface.Unfortunately, in the fra 445 of the T02/TS05 approach it would require one to iteratively recalculate the entire set of def parameters, a computationally unfeasible task.Another deficiency of the deformation meth it works fairly well only for the lowest 1st and 2nd Fourier harmonics, but rapidly deterio higher-order terms, which does not allow one to model the magnetic effects of azimuthally Birkeland currents.Also, it is a priori assumed in the model that the currents have no a 450 component, which prevents one to explore the FAC closure in the distant magnetosphere.
An alternative way to build a model with much greater flexibility is to evaluate the magn B FAC due to FACs using Biot-Savart integration.In terms of a vector potential, the problem to calculating  approach was developed independently by Ontiveros et al. (2006) in their model of the R2 FACs and the magnetopause currents.As an illustration, Fig. 6 presents a sample test distribution of the model R1 and R2 FAC density at the ionospheric level, obtained using the above summation procedure.
A great advantage of the Biot-Savart summation is that, unlike in the magnetic field deformation 470 method, it allows one to modify the current system geometry without adding unwanted artificial currents.On the negative side, in most cases the procedure is computationally rather intensive, even in its finite-segment version.

Geodipole tilt effects.
The Earth's dipole axis is inclined by ≈ 10 • (as of the 2010 epoch) to its rotation axis, which,  all current systems and lead to major deformations of the entire magnetospheric configuration.In the inner magnetosphere, the spatial distribution of trapped particles is controlled by the strong internal field, so that the ring current is nearly rigidly "attached" to the solar-magnetic (SM) equatorial plane.At distances larger than the "hinging radius" R H ∼ 8 R E , the effect of the solar wind entrainment comes into play, which makes the equatorial current sheet gradually deflect away from the SM equatorial plane.In the distant tail the current sheet aligns parallel to the solar wind flow, while at intermediate distances it bends in the form of a troughed surface.The tail current sheet deformation was thoroughly studied in the past, starting with the early work of Russell and Brody (1967).A detailed list of references can be found in (Tsyganenko and Fairfield, 2004).
A similar effect was recently discovered in the magnetosphere of Saturn, where the tilt-related deformation was shown to exist not only in the magnetotail, but also on the day side (Arridge et al., 2008).Owing to fast rotation of the planet, the Kronian equatorial current sheet takes the form of a thin disk, extending over the entire 360-degree range of longitudes.It was found that the planetary dipole tilt results in a bowl-shaped deformation of the equatorial current disk, with its periphery deflected in the direction of the solar wind velocity component, normal to the dipole equatorial plane.At an intuitive level, this can be likened to a kind of "blowing away" of the distant current sheet by the incoming solar wind, even though the solar wind does not actually penetrate into the magnetosphere.Although the relative magnitudes of electrodynamic and centrifugal forces in the Earth's case are quite different from those at Saturn, the basic physics of the tilt-related warping of the equatorial current should be the same.Therefore, this finding may provide helpful insights into the mechanism of current sheet deformation or, at least, suggest an optimal mathematical form of its empirical description.
As already noted above, however, the hardest problem in the empirical modelling is not to describe the electric current geometry, but to represent the associated magnetic field.In the case of dipole tilt effects, it is natural to begin with an untilted configuration with = 0 and use it as a starting approximation; the central question here is how to extend it to the tilted case.An effective answer and a powerful tool is the field deformation method (Stern, 1987;Tsyganenko, 1998), whose essence consists in a suitable modification of coordinates entering as arguments in the original vector field, transforming the latter into the desired final configuration.
The approach adopted in the most recent models is to start with untilted, fully shielded symmetric configurations and then apply two consecutive deformations.The first one is a rotational deformation around the Sun-Earth axis, resulting in a trough-like warping of the current sheet in the Y − Z plane.Owing to the axial symmetry of the undeformed magnetopause, its original shape remains intact at this step.The second deformation is a radially dependent rotation around the Y GSM axis by the angle * (r), so that * (r) ≈ in the inner magnetosphere, but gradually falls off to zero in the distant tail.As a result, the equatorial current sheet follows the dipole equatorial (solar-magnetic) plane in the inner magnetosphere, then bends at r ∼ R H and gradually becomes parallel to the GSM equatorial plane at Z GSM ≈ R H sin , in agreement with observations.Mathematically, in both cases the corresponding modification of the magnetic field is accomplished by a two-step procedure.First, the old (undeformed) coordinates r are replaced by new ones r = r (r) in the original equations for the undeformed field components, which yields "interim" field components B * (r) = B[r (r)].Second, final deformed field components are obtained as B = T • B * , where the tensor T is composed of partial derivatives of the components of the r vector with respect to those of r.Details of the method can be found in the original papers cited above.
In recent empirical models, the deformation procedure described above was applied to untilted shielded fields of all sources residing inside the model magnetopause, including both the equatorial and field-aligned currents.However, there is a subtlety.Since the rotational "spacewarping" is applied to the entire shielded field, it deforms not only the currents inside the magnetosphere, but also the shape of the model magnetopause, and the problem is to make the deformation consistent with independent data on the position of both the equatorial current and the magnetospheric boundary.To partially mitigate that problem, in the T02, TS05, and TS07D models the rotational deformation was modified in order to bring the magnetopause tilt-related shift in agreement with its statistically observed amplitude in the tail.However, that method is not yet flexible enough to reconcile the bowlshaped deformation of the entire equatorial current with the global tilt-related deformation of the magnetopause, built in its most recent models, such as that by Lin et al. (2010).In the T13 model, an attempt has been made to employ a more sophisticated and flexible tilt-related spacewarping, optimized by means of a joint fitting procedure, so that both the magnetopause and the equatorial current sheet deform in the required manner.

Magnetopause currents and the shielding field
Owing to its curl-free nature inside the magnetosphere, the magnetopause field is commonly represented as the gradient of a scalar potential B MP = −∇U , satisfying Laplace's equation ∇ 2 U = 0 with the Neumann boundary condition , where B i is the field of an intramagnetospheric source to be shielded.As already noted in Sect.3, the total shielding field is usually split into a linear combination of partial fields h i , corresponding to each term in the expansion in Eq. ( 6).Specific forms of the scalar potentials u (k) i depend on the geometry of a source to be shielded, and the most effective and commonly used method to derive their parameters is based on minimizing the residual rms value of the normal component over the part of the boundary confining the modelling region (Tsyganenko, 1995).In its original form, the method dates back to the "source-surface" approach by Schulz and McNab (1987).
The shape and size of the model magnetopause are a priori described by an analytical surface, whose parameters are functions of the solar wind ram pressure, and, in more recent models (e.g., Shue et al., 1998), of the B z -component of the IMF.The most sophisticated recent model by Lin et al. (2010) is also parameterized by the dipole tilt angle.
There exists a great variety of methods to compose suitable shielding fields.An abundant source of scalar potentials is a suite of solutions of Laplace's equation in several coordinate systems that allow separation of variables (e.g., Moon and Spencer, 1971).A general approach here is to choose solutions taking into account the geometry of a source to be shielded, its parity, and asymptotic properties.For example, it is a priori clear that the spherical harmonic expansions in negative powers of r, Eq. ( 2), used for the IGRF model of the Earth's main field, is a poor choice as a shielding field for the external sources, because they diverge at r → 0. Likewise, similar expansions with positive powers of r are equally unsuitable because of their divergence at r → ∞.By contrast, the scalar potentials obtained in parabolic (Alexeev and Shabansky, 1972;Stern, 1985) or cylindrical (Beard et al., 1982;Tsyganenko, 1995) coordinates, are a much better option, owing to their gradual monotonic variation along the Sun-Earth axis.
In many cases the residual rms σ i on the boundary can be significantly reduced by including the fields of static image sources (such as straight line currents, loops, current sheets, etc.), placed outside the modelling region and optimized by varying their magnitudes and geometrical parameters.An old archetype example is the image dipole model of the magnetopause field (e.g., Taylor and Hones, 1965;Antonova and Shabansky, 1968), providing a fairly accurate and extremely simple approximation of B MP in the inner magnetosphere.However, in modern global empirical models the fields of image sources are used only as supplementary terms, added to the main expansions to further improve the shielding quality.
The most recent models (T02, TS05, TS07) have employed probably the simplest and most effective version of the shielding field, based on solutions of Laplace's equation in Cartesian coordinates.These so-called "box" (or rectangular) potentials, first used in the T96 model for shielding the field of the tail current sheet, exponentially decrease tailward and contain sines and cosines of the scaled coordinates Y and Z.A commonly used version of the shielding potential is a linear combination of the box harmonics with N 2 coefficients a ik and N nonlinear parameters p i .The expansion in Eq. ( 26) yields a magnetic field whose symmetry properties correspond to the case of an untilted geodipole, that is, both B x and B y components are odd with respect to Z, while B z is even.As discussed above in Sect.3.3, the tilt-related deformation of the field of the equatorial and field-aligned currents is modeled by applying a spacewarping procedure to the untilted symmetric shielded field.It is natural therefore to use Eq. ( 26) as a universal generic form of the shielding field for those sources.
The situation is somewhat different in the case of shielding the Earth's main field, because the magnetic moment of a tilted dipole can be split into a sum of components, parallel and perpendicular to X GSM axis This allows one to conveniently represent the scalar potential of the total shielding field as a sum of two independent solutions, weighted by sin and cos where U ⊥ is represented by Eq. ( 26) and The above expansion for U is similar to Eq. ( 26), except for its parity with respect to Z.Here B x and B y are even, and B z is an odd function of Z.This allows us to avoid using the deformation method and to obtain the dipole shielding field with great accuracy by simultaneously deriving the optimal coefficients a ik , b ik , and the nonlinear parameters p i , q i in a least-squares minimization procedure.Figure 7 shows sample configurations of the shielded field of Earth's dipole, for its untilted and tilted orientation.
A separate important question is how to take into account in the shielding field the effects of variations of the magnetopause size and shape, induced by changing solar wind conditions.A common solution is to re-scale the shielding fields and/or represent the coefficients and nonlinear parameters entering in Eqs. ( 26) and (29) as polynomials of the corresponding driving parameters.More details on that are given in Sect.4.2 below.

Parameterization of the empirical models
Fitting the empirical expansions in Eq. ( 6) to the entire body of spacecraft data would provide only an average model configuration, without any information on the response of the magnetosphere to changing interplanetary conditions.Models of that sort are, however, of little value, since the main goal of the modelling is to reproduce the dynamics of stormtime space weather events.The hardest problem is the extreme disparity between the enormous multitude of possible disturbance scenarios and the fact that, at any specific moment in time, the magnetosphere is monitored by no more than just a few spacecraft.In most cases, their observations are supported by the simultaneous data of solar wind probes and ground-based geomagnetic observatories (e.g., in the form of activity indices).In some cases, low-altitude data on the particle precipitation boundaries are also available.However, these sparse data are obviously insufficient to faithfully reproduce the instantaneous geomagnetic field structure.The nents, parallel and perpendicular to X GSM axis This allows one to conveniently represent the scalar potential of the total shielding field as a sum of two independent solutions, weighted by sin Ψ and cos Ψ 590 where U ⊥ is represented by Eq.( 26) and The above expansion for U is similar to Eq.( 26), except for its parity with respect to Z.Here B x and B y are even, and B z is an odd function of Z.This allows us to avoid using the deformation method 595 and to obtain the dipole shielding field with great accuracy by simultaneously deriving the optimal coefficients a ik , b ik , and the nonlinear parameters p i , q i in a least-squares minimization procedure.
Figure 7 shows sample configurations of the shielded field of Earth's dipole, for its untilted and tilted orientation.
A separate important question is how to take into account in the shielding field the effects of vari-600 ations of the magnetopause size and shape, induced by changing solar wind conditions.A common solution is to re-scale the shielding fields and/or represent the coefficients and nonlinear parameters entering in Eqs.( 26) and ( 29 principal goal and central idea of the data-based modelling is to engage vast information contained in the archived data, to effectively improve and maximize the accuracy of the magnetic field reconstruction in specific events of interest. Two different approaches can be envisaged, the first of which is the data sorting method.In essence, it is based on establishing quantitative criteria and choosing appropriate parameters to uniquely select from the entire grand data archive smaller subsets, obtained under conditions similar to those that existed during the modeled event.Model configurations derived from the subsets reflect to some extent the average trends in the magnetic field restructuring in response to changes in the controlling parameters.By a proper choice of the selection criteria, the size of individual subsets can be optimized to reach a trade-off between the accuracy and the resolution of the modelling.In a primitive form, that approach was used in early models, based on binning of then available data into intervals of the Kp-index (Tsyganenko, 1990, and references therein).Recently, an advanced dynamical data selection method was developed (Tsyganenko and Sitnov, 2007;Sitnov et al., 2008), in which not only a groundbased index (SYM-H), but also the solar wind/IMF data are taken into account.This so-called "nearest-neighbour" approach was realized in the TS07D model, successfully used in the empirical reconstruction of the storm-time evolution of the magnetosphere during specific events (Sitnov et al., 2010(Sitnov et al., , 2012)).
The second approach seeks from the outset to relate the magnitude and geometry of individual field sources with geoeffective characteristics of the incoming solar wind and ground activity indices by means of "quasi-universal" equations, whose a priori unknown parameters are derived once and for all from the entire grand set of archived data.Historically, the first model of that kind was T96, in which the magnitude coefficients of the magnetopause, ring, tail, and field-aligned current modules were represented as functions of hourly averages of the solar wind pressure, IMF, and the Dst-index.However, the model did not take into account previous solar wind conditions, nor any effects of magnetospheric inertia.More recent models were based on data with finer (5 min) resolution, and first attempts were made to reproduce the delayed reaction of the magnetosphere.Thus, in the T02 model the external driving parameters were averaged over one-hour intervals preceding the current time moment.A more advanced method was implemented in the TS05 model, in which the magnitude coefficients were represented as solutions of equations, empirically approximating the dynamics of individual field sources.In following Sects.4.1 and 4.2 we discuss both approaches in more detail.

Binning by the Kp-index in early models
Due to lack of continuous solar wind and IMF observations in the beginning of space era, the only available data that could be used to quantify the state of the magnetosphere came at that time in the form of ground-based activity indices.Accordingly, the first empirical model based on space magnetometer data (Mead and Fairfield, 1975; MF75 for short) as well as the TU82, T87, and T89 models, were parameterized by the Kp-index.A standard approach was to group all data into several subsets, corresponding to consecutive bins of Kp, and to separately derive model parameters for each subset.The model configurations revealed a systematic increase of the external magnetic field and associated currents with growing Kp, as illustrated in Fig. 8 based on T89 model.The prominent peaks of the azimuthal current volume density around midnight are due to the fact that the equatorial current sheet is thinnest at Y = 0 and expands towards the tail flanks.Due to the crude nature of the Kp-index, derived from three-hour fluctuation amplitudes of the ground magnetic field (Bartels et al., 1939), the above models could not replicate the actual dynamics of the magnetosphere, a very complex system with finite response/relaxation and loading/unloading timescales.It should also be kept in mind that even more advanced indices provide only "integrated" information on the external currents, in which contributions from all individual sources are mixed up.Finally, as already noted in Sect.2, an inevitable restraint inherent to all empirical models is that they are based on sets of asynchronous observations made at a vast range of locations at different times and during events with largely different time histories.

"Nearest-neighbour" method and the TS07D model
The central idea in this approach is to generalize and refine the binning method by introducing a more sophisticated set of variables, quantifying the evolution of the magnetospheric state during an event, and use it to define a "similarity" criterion for selecting data records from the grand data archive (Sitnov et al., 2008).In a more rigorous formulation, the criterion reduces to the requirement that the normalized "state vector" G N N , corresponding to the i-th "nearest neighbour" data point, falls into a limited neighbourhood of the state vector G, corresponding to the current data record, as illustrated in Fig. 9.In the TS07D model, the state vector G = { vB z , SymH , DSymH/Dt } had three components, defined as 6-h averages of the solar wind electric field, SymH-index, and its time derivative, respectively.Such a choice was motivated by the fact that the above parameters are principal variables, defining (in terms of the Dst-index) the state of the magnetosphere in the well-known equation by Burton et al. (1975).The averages were centred at the current time moment and normalized by standard deviations of the corresponding quantities.
As the state vector G and the binning "sphere" move in the parametric space with time, some "neighbour" points (data records) exit from the subset, while new points enter in.As a result, the sliding selection procedure generates a sequence of subsets, covering the time interval of interest.Fitting the model to each subset in the sequence yields consecutive sets of the model parameters and field configurations, representing the dynamics of the magnetosphere during the event.Figure 10 shows three distributions of the equatorial electric currents, derived from the TS07D model field for the storm of 21-23 April, 2001.The plots correspond to the early main phase (panel a), the peak of the SYM-H index (at ∼ −100 nT), and the late recovery phase (panel c).The modelling reveals the formation of a hook-shaped partial ring current during the main phase and its decay during the recovery phase, with the formation of an extended axisymmetric ring current.This demonstrates the great potential of the approach, especially in view of the continuing rapid inflow of available data from ongoing and future magnetospheric single-and multi-spacecraft missions.

Parameterization by "global" driving variables
By the beginning of the 1990s, a sufficiently large amount of archived interplanetary data was accumulated and made available (King, 1994).Fairfield et al. (1994) compiled a large set of magnetospheric magnetic field data for the period 1966-1986, tagged by hourly averages of the solar wind plasma and IMF data.The data set was used to calibrate the T96 model and the OM97 model by Ostapenko and Maltsev (1997), which differed from each other in two aspects.First, the OM97 model (like the MF75) described the net external field vector by a single set of formal expansions in powers of coordinates and driving variables, and, Fig. 9. Illustrates the "nearest-neighbour" selection of archived data into a binning sphere, drifting in the parametric space (Sitnov et al., 2008).
for that reason, its region of validity was limited to the inner magnetosphere.Second, like all other early models, it had no explicit magnetopause.By contrast, the T96 model was intended as a global model and employed the modular approach, in which all individual sources were separately shielded inside a predefined magnetopause, with their magnitude being driven by the hourly Dst index and concurrent interplanetary parameters.The T96 ring current was parameterized by representing its magnitude coefficient a RC as a linear function a RC = a RC,0 + a RC,1 Dst * (30) of the "corrected" Dst-index Dst * = 0.8 Dst − 13 P dyn , where the coefficient 0.8 compensates for the amplification of the H-component of the ground disturbance field due to the induction currents inside Earth, and the second term removes the variable contribution from the magnetopause currents.The assumed value 13 of the coefficient of P dyn was derived from a pressure balance equation for the adopted shape of the model magnetopause (Tsyganenko, 1996).Note, however, that in a recent study by Zhao et al. (2011), a strong dependence of that coefficient on the disturbance intensity was found.
The magnitude coefficients of the modules representing contributions from the tail current were assumed to have the generic form where the driving parameter in the last term γ = P 1/2 dyn B t sin(θ/2) included the transverse component of the IMF B t = (B 2 y + B 2 z ) 1/2 and its clock angle θ.The assumed form in Eq. ( 32) quantifies the fact that the solar wind ram pressure and the reconnection with the IMF are principal factors defining the magnitude of these sources.A similar form (but with a 1 = 0) was also adopted for the FAC modules.
A separate important question is how to model the effects of the solar wind pressure and IMF variations in the shielding component h (k) i , entering in each module according to Eq. ( 6).Regarding the pressure effects, the task is somewhat facilitated by the fact that, according to recent magnetopause models (Shue et al., 1998;Lin et al., 2010), the average boundary responds to the pressure variations by expanding/contracting in a self-similar way, i.e., without changing its shape.The calculation of modified field vectors becomes especially simple for the dipole shielding field, B MP,dip , since the dipole field itself is self-similar, which results in a compact re-scaling equation from a standard pressure P to its new value P as B MP,dip (r, P ) = κ 3 B MP,dip (κr, P ), where κ = (P /P ) α and the power index α is close to its theoretical value 1/6.Note that, according to Lin et al. (2010), α ≈ 0.19.The case of non-dipolar sources (RC, TC, and FACs) is more complex.In general, there is no reason to assume that their average geometry re-scales self-similarly in response to the magnetopause compression/expansion.In the T96 model, nevertheless, self-similarity of all current systems was assumed from the outset, to avoid complications with the shielding field derivation.All later models also retained that assumption, except the most recent T13, in which individual modules are independently scaled, regardless of the magnetopause size.
Concerning the dependence of the magnetopause shape on the IMF and its impact on the shielding field, the principal fact is that, according to existing boundary models by z ] 1/2 , normalized by its rms magnitude (∼ 5 nT).Based on this, an extended set of coefficients is obtained by minimizing σ i in Eq. ( 25) over a cumulative set of boundary locations, generated for several values of IMF B z , evenly distributed within a reasonable range.That approach was implemented in the T13 model and yielded satisfactory results.

Dynamically driven models
A fundamental paradigm at the core of data-based modelling of the magnetosphere consists in two assumptions.First, the wide range of possible magnetospheric configurations can be defined with sufficient accuracy by specifying magnitudes and geometrical parameters of a limited set of modules, representing contributions to the total field from individual current systems.Second, there exists a deterministic relationship between the above parameters and routinely monitored internal and external observables, quantifying the current state of the magnetosphere and the interplanetary medium, as well as their previous history over an interval of up to a few days.The second assumption is equivalent to the statement that global field configurations, corresponding to different events but with similar patterns of external input and internal conditions, are also similar to each other.In fact, the above postulates serve as the basic justification for using asynchronous data archives for reconstructing the dynamics of specific events.It should be understood however, that the assumed determinism is naturally limited by inherent chaotic processes and instabilities, resulting in a dramatically large variety of actual magnetospheric structures and their evolution scenarios, which are extremely hard to predict even with the fullest knowledge of the solar wind and IMF conditions.
The essence of the "dynamically driven" approach is to combine an empirical model of the spatial distribution of the magnetic field with empirical equations, which relate the temporal behaviour of individual current systems to the external driving.The treatment is based on an assumption that each magnetospheric current system has two types of response to the external driving.The first is associated with disturbances due to variations of the solar wind pressure, rapidly propagating via Alfvén waves inside the magnetosphere.On the timescale of storms, this is a virtually instantaneous reaction, which can be easily reproduced by including an appropriate pressure-dependent factor in the size of the magnetopause and the related strength of the Chapman-Ferraro field.The second type of response is associated with slower processes, such as reconnection at the magnetopause, plasma convection, particle losses due to pitch angle diffusion, charge exchange, etc.These effects can be empirically modeled by including in the magnitude coefficient a (k) i in Eq. ( 6) a term W , varying in time according to the equation Here the first "source" term S(t) on the right-hand side represents the external driving.Its specific form can be assumed similar to one of the solar wind-magnetosphere coupling functions (e.g., Newell et al., 2007, and references therein).In the TS05 model, the source term was assumed to be a product of powers of three principal interplanetary medium parameters, affecting the magnetospheric state, S = aN λ V β B γ s , where N, V , and B s are the solar wind density, speed, and the southward IMF component, respectively, raised to a priori unknown powers λ, β, and γ , to be derived from data.
The second "loss" term L(W ) represents the decay rate of the field source.Its physical interpretation depends on which current system is being considered.For example, in the case of the ring current it is related to the dissipation of the energetic particle population due to charge exchange, precipitation, drift losses, etc.In the TS05 model, it was assumed to be proportional to the excess of W over its quiet-time value W 0 : L(W ) = r(W − W 0 ).In this case, Eq. ( 33) becomes similar to the well-known linear equation of Burton et al. (1975), which yields an exponential relaxation of W to its pre-storm level after the external driving is turned off.It is worth noting that, according to Aguado et al. (2010), the magnetosphere recovers faster during first hours of the recovery phase than at later times, which implies a hyperbolic relaxation rate with L ∼ W 2 .Figure 11 (from Tsyganenko and Sitnov, 2005), shows in the top panel the dynamics of the observed and predicted variation of the SYM-H index (in the original paper the vertical axis was erroneously labeled as Dst), inferred from the TS05 model for a 12-day interval of a long double storm of 3-14 September, 2002.The overall agreement between the observed SYM-H (heavy black line) and that derived from the model (thin black line) is remarkably good in this case.The six coloured traces correspond to separate contributions to the index from individual current systems, as explained by the legend.In agreement with a conjecture by Alexeev et al. (1996), the principal contribution to the SYM-H/Dst index at the peak of the storm main phase comes not from the SRC alone, as was believed to be the case since the early days of space physics, neither solely from the TC, as was argued by Maltsev (2004), but in roughly equal  (Tsyganenko and Sitnov, 2005).
were found to differ significantly from each other, from as large as ∼30 hours for the SRC to only ∼50 min for the R1 FACs.The total magnitude of the currents were also found to vary dramatically in the course of major events, with peak values as large as 5-8 MA for the SRC and R1 FACs.In line  (Tsyganenko and Sitnov, 2005).
shares from both these sources.However, due to its much shorter relaxation time, the TC contribution decreases much faster than that of the SRC during the recovery phase.
Another interesting piece of information provided by the TS05 model concerns the dynamics and peak values of the total current corresponding to individual field sources.Their relaxation/response timescales were found to differ significantly from each other, from as large as ∼ 30 h for the SRC to only ∼50 min for the R1 FACs.The total magnitude of the currents were also found to vary dramatically in the course of major events, with peak values as large as 5-8 MA for the SRC and R1 FACs.In line with simulation results by Liemohn et al. (2001), at the peak of the main phase the total PRC can largely exceed the SRC, reaching 10 MA and even more, but it quickly subsides, as the external solar wind driving disappears, with the relaxation time less than 2 h.The TC increases dramatically during the main phase and shifts earthward, so that the peak current concentrates at unusually close distances of 4-6 R E .

Spacecraft data for modelling
Space magnetometer data, complemented with concurrent interplanetary and ground-based observations, are one of the cornerstones of empirical modelling.Since the beginning of the space era, huge archives have been created, containing an enormous body of data taken by many satellites at different locations, seasons, solar cycle phases, and disturbance levels.Mead and Fairfield (1975) compiled the first set of distant magnetospheric magnetic field data, taken by four IMP spacecraft during 1966-1972 and used it to create the MF75 model.Tsyganenko and Usmanov (1982) added HEOS-1 and -2 data to the Mead-Fairfield set and developed a more realistic TU82 model with explicitly defined ring and tail current sources.The data set was further extended by Tsyganenko and Malkov (described by Peredo et al., 1993), who added ISEE-1 and -2 data from 1977-1981, while Fairfield independently added HEOS observations and additional IMP-6 data to the original Mead-Fairfield database.Editing those data and merging them into one large database resulted in a data set covering the period 1966-1986, described by Fairfield et al. (1994) and subsequently used in the derivation of the T96 magnetospheric magnetic field model.
In the following years, the launch and prolonged operation of AMPTE/CCE/IRM (1984-1988), Geotail (1992-present), Wind (1994-present), Polar (1996Polar ( -2008)), ACE (1997present), Cluster (2001-present), Themis (2007-present), as well as the succession of geosynchronous GOES satellites, resulted in rapid, manifold, and continuing expansion of the available database.It is hard to overstate the impact and challenge of that wealth.The abundance of the data has been instrumental in the development of empirical field models, since the accuracy of the latter critically depends on the range and density of data coverage not only in geometric, but also in parametric space, including the dipole tilt angle and solar wind/IMF parameters.
It should be realized that most of the available data correspond to quiet and weakly disturbed time intervals, while unusual and strongly disturbed periods (most interesting for the physics and most important for space weather) are relatively rare, so that storm-time data constitute only a few percent of the entire database.The ultimate goal of the modelling is to fairly reproduce the entire range of magnetospheric states, covering both undisturbed periods and all disturbance phases.Hence, the modelling data sets should be compiled in such a way that they contain nearly balanced amounts of quiet pre-storm and storm-time data.An optimal method to construct a data set is to organize it as a collection of events, each of which starts with a one-or two-day period of quiescence, followed by a disturbed period of nearly the same or longer duration (Tsyganenko et al., 2003).
Another important issue is the choice of time resolution of the data in the modelling sets.A 5 min time interval corresponds to a ∼ 20 R E travel distance in the solar wind flow, which is commensurate with the transverse scale size of the magnetosphere.It therefore can be adopted as a characteristic timescale for the magnetosphere to respond to changes in the external pressure.Adopting finer resolution would lead to unreasonably large and redundant data sets, while longer averaging intervals would smear out short-term variations of the field due to incoming shock fronts or transient gusts of the solar wind.Also, in view of the large distance between the location of the interplanetary medium monitors (usually, at the L1 point with X GSM ∼ 220R E ) and the subsolar bow shock, it is hard to evaluate the solar wind travel time with an accuracy significantly better than a few minutes, even using sophisticated propagation techniques, e.g., like that employed by Weimer et al. (2008).One more argument in favour of the 5 min average data is that concurrent interplanetary parameters are routinely available from the OMNI website with this time resolution.
Finally, the events included in the modelling sets should be continuously covered by concurrent interplanetary medium data.Short (less than a few hours) gaps in those data can be filled by interpolation, but only if there is no indication of major changes in the solar wind, which in most cases is evident from the absence of irregularities in the SYM-H index.Data taken during disturbed intervals, but not covered with simultaneous interplanetary observations are absolutely useless, so it makes no sense to include such events in the modelling set.
Routine procedures typically involved in the data processing include initial retrieval of fine-resolution data, their reformatting and removal of the internal (IGRF) part from the total field vectors.Then the high-resolution data are averaged over 1 min intervals and merged with concurrent interplanetary data, which allows us to make an initial automatic removal of data taken outside the magnetosphere, using a magnetopause model driven by solar wind parameters, e.g., that by Lin et al. (2010).The 1 min average data are then visually inspected day-by-day, in order to eliminate bad/questionable data records and those clearly belonging to remaining unfiltered magnetosheath intervals.In addition, all data taken at geocentric distances r ≤ R c are removed, where the inner limit R c is typically between 2.5 and 3.5R E .Such nearperigee data are usually contaminated by large errors, rapidly increasing with decreasing r due to the fast-growing main geomagnetic field and inaccuracies of the satellite attitude data.The final step is to average the external field components over 5 min intervals.
To illustrate the overall coverage of the magnetosphere by space magnetometer data, Fig. 12 displays the spatial distribution of observations in a grand modelling set containing 241 138 data records, used in the calibration of the most recent T13 model.The set includes 123 storm events during the period 1997-2012, with peak negative values of SYM-H not exceeding −200 nT.
A separate problem is a strongly nonuniform spatial distribution of the data, with much fewer data in the middle and distant tail due to much lower number of high-apogee spacecraft and their much longer orbital period, in comparison with the inner magnetospheric missions.Figure 13 displays in the upper panel the radial distribution of data in the grand set shown in Fig. 12 (note the logarithmic scale of the vertical axis).It is clearly seen that the largest portion of the data is confined within r ≤ 12 R E , owing to the strong disparity between the relatively sparse population of Geotail and Cluster data points in the midtail region and a much denser coverage of the inner magnetosphere by Polar and Themis.Another data disparity factor, illustrated in the centre panel of Fig. 13, is the rather steep decrease of the magnetic field magnitude from the inner magnetosphere to the distant tail.
In this situation, using unweighted data in the model calibration might result in a significant bias of the reconstructed field in the underpopulated magnetotail, where, in addition, the field itself is relatively weak and, hence, sensitive to even small changes of best-fit values of the model parameters.This can be avoided by introducing a weighting procedure, in which each consecutive adjacent r = 0.5 R E bin of radial distance is given a weight W i = ( N /N i )( B / B i ), inversely proportional to the number of data records and to the average field magnitude in that bin.The bottom panel in Fig. 13 shows the radial distribution of the weight function, N. A. Tsyganenko: Data-based models of the magnetosphere is evident from the absence of irregularities in the SYM-H index.Data taken during disturbed intervals, but not covered with simultaneous interplanetary observations are absolutely useless, so it makes no sense to include such events in the modelling set.Routine procedures typically involved in the data processing include initial retrieval of fine-resolution data, their reformatting and removal 860 of the internal (IGRF) part from the total field vectors.Then the high-resolution data are averaged over 1-min intervals and merged with concurrent interplanetary data, which allows us to make an initial automatic removal of data taken outside the magnetosphere, using a magnetopause model driven by solar wind parameters, e.g., that by Lin et al. (2010).The 1-minute average data are then visually inspected day-by-day, in order to eliminate bad/questionable data records and those clearly 865 belonging to remaining unfiltered magnetosheath intervals.In addition, all data taken at geocentric distances r ≤ R c are removed, where the inner limit R c is typically between 2.5 and 3.5R E .Such  A separate problem is a strongly nonuniform spatial distribution of the data, with much fewer data in the middle and distant tail due to much lower number of high-apogee spacecraft and their much longer orbital period, in comparison with the inner magnetospheric missions.Figure 13

Assessing the performance of empirical models
When comparing the models with each other and evaluating their accuracy, one should keep in mind the inevitable mismatch between, on the one hand, the complexity of the real magnetosphere, its vast dimensions, the broad range of spatial/temporal scales involved, and the very wide variety of possible event scenarios, and, on the other hand, the inherent limitations of the model description of the field structure and dynamics.For that reason, it is virtually impossible to create a "universal" global model, equally accurate at any distance and for any conditions.Also, the very notion of accuracy may imply different meanings, depending on the model's specific application.For example, a model can be used either for mapping the geomagnetic field lines between a spacecraft location and the ionosphere (that issue is discussed below in Sect.6.1), or to trace energetic particle orbits in the geomagnetic field.In the latter case, the principal quantity of interest is the magnetic field vector B, and the overall accuracy of the model can be best estimated (and optimized) in terms of the rms deviation of the model field from data, normalized by the rms magnitude |B| of the observed field.Another possible way to assess a model's quality is to calculate the correlation coefficients between the observed and model field components.As an example, Fig. 14 shows scatter plots of the observed vs. model field components for the most recent T13 model.Three plots at the top of the figure correspond to the total field, including the Earth's contribution.Here all points with r ≤ 6.6 R E were excluded from the comparison; due to the overwhelming dominance of the internal field in the inner magnetosphere, adding those data would further improve the correlation, but would suppress and hide information on the external field model performance.To better highlight the latter, the three plots at the bottom show the same data points, but only for the external field components, i.e., with the IGRF part subtracted.
When estimating the agreement between the model and observed fields, it is more convenient to use a single vector correlation coefficient R v , instead of the three separate coefficients for B x , B y , and B z .The R v coefficient is defined only by the mutual orientation of the corresponding individual vectors in the data set and, hence, is independent of the choice of the coordinate system.Formally, it has the same properties and is defined in exactly the same way as the commonly used correlation coefficient for scalar data, i.e.
except that the scalar quantities are replaced here by a set of vectors B (obs) i and B (mod) i , representing the observed and model external fields, respectively.
In the case of field line mapping, the most relevant quantity is the total magnetic field direction , representing the observed and model external fields, respectively.
In the case of field line mapping, the most relevant quantity is the total magnetic field direction vector b = B/B, entering in the field line equation dr/ds = b(s).This suggests deriving model parameters by minimizing the rms deviation between the observed and model b vectors, instead of that for the full B vectors.That approach was implemented in the derivation of the T96 model parameters and yielded quite robust tail field configurations.However, at closer geocentric distances the merit function based on the directional criterion becomes progressively less and less sensitive to the external field, because of the rapid growth of the Earth's main field.In the inner magnetosphere the geomagnetic field remains nearly quasi-dipolar at almost all times, except during strong storms.That could be the most likely cause of the overstretched T96 model field in the inner magnetosphere (Tsyganenko, 2002b;McCollough et al., 2008).
The quality of a model from the mapping viewpoint can be quantified by calculating histograms of the angular difference between the observed and model B vectors over a sequence of bins of radial distance.An example is shown in Fig. 15, where the left and right panels correspond to the inner and outer magnetosphere, respectively.Each plot displays three histograms, for the T96 (red) and T13 (blue) models, and for the internal field (IGRF) without any external field model (black).As can be seen from the plots, the most recent model yields the best result, though with only a token improvement in the distant magnetosphere.

Evaluating the mapping errors
The histograms in Fig. 15 provide distributions of the local angular deviations between the observed and model field vectors.However, the quantity of ultimate interest in most mapping applications is the integral error of a field line footpoint location.A method to estimate those errors was suggested by Pulkkinen and Tsyganenko (1996), based on a first-order perturbation technique.Its essence is to evaluate the cumulative equatorial shift R eq of a model field line L m from the "actual" one L a that starts from the same ionospheric footpoint.Assuming that both field lines are not too far from each other, the net shift can be found by integrating local deviations δs = δbds along the model field line between the ionosphere (s = 0) and the equatorial plane (s = S eq ), so that  defined using data points inside a field line tube of a finite thickness, centered on the field line L m .The above method was used to test the T89 model, with results presented in the form of polar diagrams displaying 2-D distributions of the mapping errors for several levels of the Kp-index (see Plates 1-4 in the above-cited paper).The field line mapping direction can in principle be reversed, so that one could evaluate the field line footpoint deviations in the ionosphere, instead of at their apex points.
Another powerful technique to assess the models' performance in terms of the mapping accuracy is based on lowaltitude observations of energetic charged particles.The idea of the method reduces to the simple fact that, for each species of particles with a given rigidity mVc/q there exists a surface in the magnetosphere, consisting of closed geomagnetic field lines, which separates the regions with adiabatic and nonadiabatic regimes with respect to the particle's first invariant.The surface is called the isotropic precipitation boundary (IB), with its position defined by the locus of points lying on the surface of minimum B (maximum field line curvature) where the following condition is met (Sergeev and Tsyganenko, 1982;Sergeev et al., 1983) where ρ c and R L are the local field line curvature radius and the particle's gyroradius, respectively.Due to the strong magnetic field inside the boundary, γ > 8, so that particles in that region cross the minimum B region without violating their first invariant, so that the loss cone remains empty.By contrast, particles outside the IB randomly change their magnetic moments as they encounter the near-equatorial minimum B region with strongly curved field lines, and eventually end up in the ionospheric loss cone, as illustrated schematically in Fig. 16.The first factor in Eq. ( 37) is fully defined by the magnetic field model, while the second factor is just the inverse of the particle's rigidity, available from low-altitude observations of the IB.This suggests using observed IB locations, monitored  by low-altitude polar-orbiting satellites, to independently test the model's mapping accuracy.Moreover, since the critical value of γ in Eq. ( 37) is energy-dependent and the latitudinal positions of IBs are routinely observed on each spacecraft pass in a wide range of energies, it becomes in principle possible to adjust standard field models so that they more accurately reproduce the actual field configuration in specific events (Sergeev et al., 1993;Sergeev and Gvozdevsky, 1995).Combining the low-altitude IB observations with simultaneous data from several equatorial satellites and properly adapted magnetic field models can further improve the mapping accuracy and reveal interesting features of the magnetospheric dynamics during disturbances (e.g., Kubyshkina et al., 2009;Shevchenko et al., 2010).
7 Consistency of the empirical B field models with magnetospheric plasma pressure Magnetospheric magnetic fields and plasmas are intimately interrelated with each other via first-principle equations, which means that empirical B field models should have some degree of consistency with observed distributions of the plasma pressure.This is especially true with respect to regions where plasma and magnetic forces are of comparable magnitude, in particular, in the nightside equatorial magnetosphere.A number of statistical studies have been made to quantify the average distribution of the magnetospheric plasma and its anisotropy (e.g., Spence et al., 1989;Lui and Hamilton, 1992;Lui et al., 1994;De Michelis et al., 1999;Tsyganenko and Mukai, 2003).
The degree of consistency of empirical model field configurations with magnetospheric plasma distributions has been discussed for many years.Walker and Southwood (1982) tested several then existing models to check whether the related magnetic stress vector [∇×B]×B is nearly curl-free, as it should be in the case of isotropic plasma pressure.Spence et al. (1987) went further and derived 2-D distributions of the plasma pressure and its anisotropy from an empirical (TU82) field model by minimizing the rms of the residual total stress vector in the radial distance range between 6 and 12 R E .Using that technique, they inferred radial profiles of perpendicular, P ⊥ , and parallel, P , pressures, based on bi-Maxwellian particle distributions.Kan et al. (1992) compared the observed profile of isotropic pressure along the tail axis with that derived by integrating the equatorial magnetic tension, calculated from an empirical (T87) field model.A similar technique was used by Kubyshkina et al. (1999Kubyshkina et al. ( , 2002) ) for testing the consistency of an event-oriented model based on a standard (T96) magnetic field, with the observed plasma pressure.The issue of equilibrium in the near-tail region was also addressed by Hesse and Birn (1993) who developed a "ballistic" relaxation algorithm to iteratively establish a force-balanced configuration, using a version of the T87 model field as a starting approximation.Horton et al. (1993) analysed the T87 field and estimated the degree of pressure anisotropy in the central midnight plasma sheet, needed for the force balance.Cao and Lee (1994) derived anisotropic pressure distributions in equilibrium with the T87 and T89 models by directly solving the stress balance equation with respect to P and P .Toffoletto et al. ( 2001) used a relaxation algorithm to derive a force-balanced 3-D magnetosphere, starting from an empirical T96 field.Cheng (1995), Zaharia andCheng (2003a, b), andZaharia et al. (2004) developed a magnetostatic code to equilibrate 3-D distributions of anisotropic plasma in the near magnetosphere, covering all local times.Using that code, they tested an empirical field model (T96) for its inconsistency with isotropic pressure and also used that model to set up boundary conditions for a selfconsistent magnetic flux function.
In the most recent study of the empirical models' consistency with the static balance assumption (Tsyganenko, 2010), a systematic method was developed to derive a distribution of ambient plasma pressure with spatially varying anisotropy, most closely consistent with a given magnetic field.The approach was to keep the original model magnetic field intact and search for an optimal distribution of anisotropic plasma, by minimizing the rms difference between the plasma and magnetic stress vectors in the leftand right-hand sides of the force balance equation over a set of points within a model plasma sheet at distances 5 ≤ R ≤ 20 R E .
Although the method can in principle be extended to the general asymmetric case, in the above study the residual force disbalance was minimized over a 2-D area in the midnight meridian plane and it was assumed that the magnetic field is dawn-dusk symmetric.For that reason, the calculations were limited to only four models, TU82, T87, T89, and T96.The anisotropy ratio P /P ⊥ was found to significantly deviate from unity in the case of TU82 and T87 models, with progressively higher values at larger tailward distances.By contrast, the more recent T89 and T96 models yielded more realistic results with nearly isotropic pressure in the tail and a moderate pancake-type anisotropy in the inner magnetosphere, consistent with observations.Figure 17 compares radial pressure profiles corresponding to a quiet-time variant of the T96 model (both P and P ⊥ ) with statistical distributions of isotropic pressure by Lui and Hamilton (1992) and Tsyganenko and Mukai (2003).For comparison, a pressure profile calculated by 1-D integration of Ampère's force along the tail axis is shown by the red trace.Even though that method yields larger pressures at all distances, the difference does not exceed ∼ 30 % at R ∼ 9 R E .

Outstanding problems and challenges
The ultimate goal of empirical modelling is to overcome limitations caused by the sparsity of simultaneous observations in the magnetosphere and make it possible to faithfully reproduce the dynamics of individual events.All the efforts of recent decades, outlined in this review, have been focused on taking the fullest advantage of the fast-growing wealth of accumulated observations and finding optimal methods to "interrogate" historical data sets.The existing rich arsenal of mathematical methods and continuing rapid progress of computer systems gives us hope that further major advances are possible.In this section we discuss some unsolved problems and critical challenges that lie ahead.

Magnetopause shape and the IMF "penetration"
The first group of problems is related to the model magnetopause.All existing models of the magnetospheric boundary are average static surfaces, obtained by fitting analytic forms to data on the magnetopause location, detected by many satellites at different times with different states of the solar wind and IMF.No dynamical features of the boundary are present in these models, so that their parameters are

Outstanding problems and challenges
The ultimate goal of empirical modelling is to overcome limitations caused by the sparsity of simultaneous observations in the magnetosphere and make it possible to faithfully reproduce the dynamics of individual events.All the efforts of recent decades, outlined in this review, have been focused on taking the fullest advantage of the fast-growing wealth of accumulated observations and finding optimal methods to 'interrogate' historical data sets.The existing rich arsenal of mathematical methods and continuing rapid progress of computer systems gives us hope that further major advances are possible.In this section we discuss some unsolved problems and critical challenges that lie ahead.

Magnetopause shape and the IMF 'penetration'
The first group of problems is related to the model magnetopause.All existing models of the magnetospheric boundary are average static surfaces, obtained by fitting analytic forms to data on the 39 tacitly treated as instantaneous functions of the external input.As a result, fluctuations of the solar wind pressure and IMF B z on a timescale of only a few minutes produce simultaneous global variations of the entire model magnetopause shape and/or size, whereas the actual process is a wave-like propagation of the disturbance from the dayside to the tail (e.g., Collier et al., 1998).At least a partial solution to the problem could be to use a "damped" (or delayed) pressure d , instead of the "instantaneous" one, P dyn , as an argument in the functional form for the model magnetopause.A convenient variant is to relate d (t) with the input P dyn (t) via the equation Setting the initial pressure d at t = 0 equal to P 0 = P dyn (0) yields a simple solution where the characteristic timescale τ controls the speed of response of d to sudden changes of P dyn .Individual shielded submodules entering in the model expansion in Eq. ( 6) can be associated with different values of the timescale τ , replicating their different rates of response to variations of P dyn in the incoming solar wind.The same approach can be implemented with regard to IMF B z , entering as a parameter in the model magnetopause equations.These modifications could help avoid unrealistically fast synchronous fluctuations of the global model field in response to short-scale variations of the interplanetary parameters, which result, in particular, in undesirable instabilities in the hybrid MHD-particle simulations of the inner magnetosphere (N.Y. Buzulukova, personal communication, 2013).
There exists another problem concerning the global models of the magnetopause.The overwhelming majority of the boundary crossings in the data sets, used to generate the existing magnetopause models, are located sunward of X ∼ −10 R E .This can be clearly seen in the plots of Shue et al. (1997), Lin et al. (2010), and Wang et al. (2013), displaying spatial distributions of their data.Only a small fraction of the data cover the tailward part of the modelling region, due to the much longer orbital period of high-apogee satellites like Geotail, IMP-8, or Prognoz, and their smaller number in comparison with the lower-apogee spacecraft.This raises a question about the accuracy of the boundary model in the distant tail.In particular, it remains unclear whether the faster expansion of the model magnetopause obtained for IMF B z < 0 persists in the distant tail, or whether this is just an artifact of unwarranted extrapolation of the model into the region with poor data coverage.MHD simulations can in principle help resolve this issue.However, at present they rather add to the confusion.In particular, Lu et al. (2011) found that larger negative IMF B z causes the tailward magnetopause to shrink, contrary to the empirical models by both Shue et al. (1998) and Lin et al. (2010).Such a strange result is most likely due to the actual boundary studied in that work being not the magnetopause proper, but the fluopause (Siscoe et al., 2001), located inward from the magnetopause and gradually converging towards the plasma sheet owing to the dawn-dusk electric field caused by IMF B z < 0.
The issue of the magnetopause shape as a function of IMF B z is quite important in global empirical modelling for the following reason.In the static perfectly shielded field configurations, a faster tailward expansion of the boundary (predicted by the magnetopause models for IMF Bz < 0) results in a more rapid tailward fall-off of the northward directed shielding field, and hence yields generally smaller B z in the plasma sheet, i.e., more stretched configurations.Conversely, a slower magnetopause flaring (corresponding to IMF B z > 0 in the magnetopause models) yields larger shielding fields on the nightside, that is, more dipole-like configurations.This is illustrated in Fig. 18, showing two fully shielded T13 model configurations, which differ from each other only by the magnetopause flaring rate.The remarkable fact is that the same effect on the tail stretching can be obtained by keeping the magnetopause shape intact, but adding a uniform field B z e z of the same polarity as the IMF B z .That way of including the largest-scale IMF effect on the empirical field was implemented in several models with a fixed shape of the an earlier interval of northward field, during which the plasma sheet fills with cool and dense plasma.
The required sequence of the IMF B z polarity, i.e., first northward, then southward, is often observed 1125 during a passage of a CME with a flux-rope structure.The opposite sequence -first southward, then northward -can result in a faster relaxation of the SYM-H index to its quiet-time level due to tail lobe reconnection in the recovery phase of the storm.The role of such a 'quenching' effect should also be studied and quantified by means of an advanced modification of the driving term in Eq.( 33).

Field-aligned currents 1130
As already noted in Section 3.2, there is considerable room for further improvements in the modelling of the global magnetic effects of the Region 1 FACs.First, more flexibility is needed to infer the actual geometry of the FAC configuration in the distant magnetosphere.In the framework of the existing sheet-like model of the Region 1 FAC, defined by Eqs.( 20)-( 23), a possible option is to unfix the parameter α in Eq.( 22) and use direct Biot-Savart integration to calculate the magnetic 1135 field of the current system.This approach has been adopted in the T13 model and yielded initial encouraging results.Figure 19 shows the evolution of the total Region 1 and 2 FACs flowing into the northern ionosphere during a storm of 8-10 March, 2008, as reproduced by the T13 model.The model currents peak at times of southward excursions of the IMF B z , preceding SYM-H dips by up to a few hours.The model yields quite reasonable pre-storm and peak values of both the Region 1 1140 and 2 total current, in good agreement with early estimates (Iijima and Potemra, 1978;Bythrow and Potemra, 1983) and more recent results (e.g., Korth et al., 2008).Note that earlier empirical models e z .In all those models, the penetration factor κ was found from data to be in the range 0.2 ≤ κ ≤ 0.8, which was then interpreted as evidence for a real penetration of the IMF into the magnetosphere.However, as it follows from the above results, the effect of the IMF B z polarity on the global degree of openness of the magnetotail magnetic flux can be at least partially reproduced in a fully shielded model in terms of the variable IMF-dependent shape of the magnetopause.Of course, there is little doubt that the real magnetosphere is open, but at this time it is hardly possible to estimate from data the penetrated magnetic flux, unless and until more accurate magnetopause models become available.Hopefully, future studies based on more data from distant high-latitude magnetosphere will make it possible to resolve this issue.

Advanced effects in the solar wind driving
Two aspects of the model parameterization need further attention.The first concerns the complexity and nonlinearity of the external driving efficiency (and, hence, of the disturbance magnitude, quantified by the peak SYM-H value) with respect to three key factors: the speed of the solar wind, the magnitude of the southward IMF B z , and the duration of the geoeffective interval (e.g., Gonzalez et al., 1994, and references therein).In this regard, the source-loss model Eq. ( 33) with relatively simple driving and relaxation terms, adopted in TS05, is only the first step in that direction.The second aspect is the role of northward IMF B z in the storm dynamics, recognized as an important factor in the "pre-conditioning" of the magnetosphere during strong CME-related events.The essence of the phenomenon, first envisioned by Thomsen et al. (2003), consists in a significant enhancement of the disturbance magnitude, if the geoeffective extended interval of strong southward IMF is preceded by an earlier interval of northward field, during which the plasma sheet fills with cool and dense plasma.The required sequence of the IMF B z polarity, i.e., first northward, then southward, is often observed during a passage of a CME with a flux-rope structure.The opposite sequence -first southward, then northwardcan result in a faster relaxation of the SYM-H index to its quiet-time level due to tail lobe reconnection in the recovery phase of the storm.The role of such a "quenching" effect should also be studied and quantified by means of an advanced modification of the driving term in Eq. ( 33).

Field-aligned currents
As already noted in Sect.3.2, there is considerable room for further improvements in the modelling of the global magnetic effects of the Region 1 FACs.First, more flexibility is needed to infer the actual geometry of the FAC configuration in the distant magnetosphere.In the framework of the existing sheet-like model of the Region 1 FAC, defined by Eqs. ( 20)-( 23), a possible option is to unfix the parameter α in Eq. ( 22) and use direct Biot-Savart integration to calculate the magnetic field of the current system.This approach has been adopted in the T13 model and yielded initial encouraging results.Figure 19 shows the evolution of the total Region 1 and 2 FACs flowing into the northern ionosphere during a storm of 8-10 March, 2008, as reproduced by the T13 model.The model currents peak at times of southward excursions of the IMF B z , preceding SYM-H dips by up to a few hours.The model yields quite reasonable pre-storm and peak values of both the Region 1 and 2 total current, in good agreement with early estimates (Iijima and Potemra, 1978;Bythrow and Potemra, 1983) and more recent results (e.g., Korth et al., 2008).Note that earlier empirical models yielded significantly lower currents, especially the T96 model.Another area of potential improvement is adding more flexibility to the local time distribution of the FAC along the Region 1 zone at low altitudes.That can be done either by including more Fourier harmonics in the longitudinal variation of the current density J , or by adopting from the outset an ad hoc non-sinusoidal function, to model its localized peaks.Finally, it remains to add more degrees of freedom to allow the model R1 FACs to divert away from the meridional planes and close either in the plasma sheet boundary layer, or via the magnetosheath.A rich experimental resource that could be very helpful in testing the low-altitude FACs configuration and dynamics in future empirical models is the AMPERE experiment data (e.g., Anderson et al., 2005Anderson et al., , 2008)).

Dayside cusps and their magnetic structure
The magnetospheric polar cusps are formed due to the specific topology of a dipole field, confined (at least, partially) within a superconducting boundary.An essential element of such a vacuum magnetic configuration is a pair of null points on the magnetopause, at which the field lines diverge and converge.In the real magnetosphere, the null points evolve into longer and wider cleft-like regions, filled with relatively cold and dense plasma of magnetosheath origin.Due to the diamagnetism of the injected plasma, the vacuum field depressions localized near the boundary expand in space and deepen in magnitude, forming extended "funnels" with a weak magnetic field, as illustrated in Fig. 20.Fairfield (1991) compared dayside observations of IMP-4,5,6, -D, and HEOS-1,2 satellites with predictions of the T87 model and found broad regions of depressed field at high latitudes associated with the polar cusps.None of the currently existing empirical field models include the polar cusp depressions.Using data of Polar satellite, Tsyganenko and Russell (1999) studied the radial, longitudinal, and dipole tilt-angle dependence of the cusp depression, and suggested a method to include it in empirical field models, based on a local stretching of the colatitude angle.
with predictions of the T87 model and found broad regions of depressed field at high latitudes associated with the polar cusps.None of the currently existing empirical field models include the polar cusp depressions.Using data of Polar satellite, Tsyganenko and Russell (1999) studied the radial, longitudinal, and dipole tilt-angle dependence of the cusp depression, and suggested a method to include it in empirical field models, based on a local stretching of the colatitude angle.The second important phenomenon associated with the polar cusps is a high correlation of the B y -component inside the cusp funnels with IMF B y .As found in a statistical study of that effect based on Polar and Cluster data (Tsyganenko, 2009), the linear regression coefficient relating B cusp y with B IMF y dramatically (five-fold!) and monotonically rises with decreasing geocentric distance from 2.4 at 6 ≤ R ≤ 7 to 12.6 at 2 ≤ R ≤ 3.Such a stunning steepness and regularity perfectly agrees with the interpretation of IMF "penetration" into the cusps in terms of the Region-0 FACs, with their orientation and magnitude being controlled by the IMF B y (McDiarmid et al., 1979;Erlandson et al., 1988;see Tsyganenko, 2009, for more references.) Both the above effects are not yet included in the models.While the cusp depression is a relatively local effect, and, as such, does not dramatically affect the global field line mapping, the transverse zonal field due to the Region 0 currents may result in significant azimuthal shifts of the field line footpoints, an interesting subject for future studies.

Concluding remarks: universal/global and specialized/local models
As already pointed out, the magnetosphere is an extremely variable object, in which a great number of possible disturbance scenarios can unfold, depending on the magnitude of the external driving, the sequence of arrivals of the incoming solar wind structures, IMF variability, etc. Notwithstanding the complexity of the system, the efforts to develop its realistic data-based models demonstrated the general robustness of the approach and its ability to faithfully reproduce not only average statistical configurations, but also the main features of magnetospheric storm-time dynamics.At the same time, given the vast spatial expanse of the magnetosphere and the lack of simultaneous in situ monitoring by a network of satellites, it is hard to expect that a single universally accurate global model will be created in the foreseeable future.The most likely and most promising routes for further progress will be to develop sets of specialized or local models, e.g., valid in more limited domains such as the inner magnetosphere, or its dayside region, or the modelling of special kind of events, like CIR-or CME-generated storms.A separate interesting and yet largely untapped area is empirical data-based modelling of substorm reconfigurations of the magnetosphere, including the dynamics of the substorm current wedge.First advances in this field have already been made (Sergeev et al., 2011(Sergeev et al., , 2013)).
Modelling the internal part B I dates back to nearly 180 yr ago when C.-F. Gauss in his seminal works laid the foundation of data-based geomagnetic field studies.His harmonic expansion in spherical coordinates {r,

Fig. 1 .
Fig. 1. (Top) Unbounded currents (red) and unshielded B field lines (blue).(Bottom) Adding the shielding field results in current closure via the magnetopause (grey shading) and fully confined magnetic field.

Fig. 2 .Fig. 2 .
Fig. 2. Electric current flow lines, corresponding to the symmetric (left) and 'quadrupole' (centre) components of the model PRC, and the resultant total configuration (right), calculated from Eqs. (8)-(11).dipolar background magnetic field Bo, which eliminated the need to numerically trace the field lines in the calculation of the electric currents from Eqs.(8-11).In addition, the axial symmetry of

Fig. 3 .
Fig. 3. Meridional sections of the axisymmetric surfaces given by Eqs.(21)-(22) for two values of parameter α.The lines are labeled by values of the footpoint latitude in degrees.surface ξ(r, θ, φ) = 0 as well as the azimuthal distribution of the FAC density.Radially independent potentials with ∂χ/∂r = 0 correspond to purely poloidal currents, flowing in the meridional planes φ = const.Introducing a radial variation in χ adds an azimuthal component to j, which can be used 410

Fig. 3 .
Fig. 3. Meridional sections of the axisymmetric surfaces given by Eqs.(21)-(22) for two values of parameter α.The lines are labeled by values of the footpoint latitude in degrees.

Fig. 4 .
Fig. 4. Global geometry of the R1 FACs in the T02/TS05 models, shown in two projections.Left: view along the X-axis (from Sun); right: side view along the Y-axis.The electric current flow lines (yellow) close via the model magnetopause (blue shading).
for details), the procedure was somewhat different.Instead of specifying from the outset the azimuthally asymmetric FAC sheet, we started from the axially symmetric surface Θ(r, θ) = Θi0 with the constant ionospheric colatitude Θi0 of the R1 zone.The second potential χ was assumed to be a function only of the longitude φ in the simplest form f (φ) = sin mφ, with the goal to represent the local time distribution of the FACs by the first few Fourier harmonics.As in the case of the PRC (see Section 3.1.1above) the assumed axial symmetry of the surface and the sinusoidal variation of the FAC density greatly simplified the problem by making it possible to isolate the φ-dependence in the magnetic field components into separate factors sin mφ and cos mφ (Tsyganenko, 1993, Appendix B) which allowed us to reduce the problem to 2D and define the entire 3D field by calculating its components by Biot-Savart integration

Fig. 4 .
Fig. 4. Global geometry of the R1 FACs in the T02/TS05 models, shown in two projections.Left: view along the x axis (from Sun); right: side view along the y axis.The electric current flow lines (yellow) close via the model magnetopause (blue shading).

Fig. 6 .
Fig. 6.A sample distribution of model R1 and R2 FAC density (in arbitrary units) at the ionospheric level, obtained using the fast Biot-Savart integration.
475in its turn, is inclined by 23.4 • to the normal to the ecliptic plane.This results in diurnal and yearly variations of the angle Ψ between the geodipole axis and the terminator plane within the range −33.4 • ≤ Ψ ≤ 33.4 • .The dipole tilt variations affect all current systems and lead to major

Fig. 6 .
Fig. 6.A sample distribution of model R1 and R2 FAC density (in arbitrary units) at the ionospheric level, obtained using the fast Biot-Savart integration.

Fig. 7 .
Fig. 7. Sample configurations of a shielded purely dipolar field for the untilted (left) and tilted (right) case.

Fig. 7 .
Fig. 7. Sample configurations of a shielded purely dipolar field for the untilted (left) and tilted (right) case.
Shue et al. (1998) and Lin et al. (2010), negative IMF B z results in smaller standoff distances and larger tailward flaring rates of the magnetopause, while positive IMF B z leads to an opposite effect in the Shue et al. model, but causes no change at all according to Lin et al.These effects can be taken into account by representing the coefficients a ik in Eq. (26) as Taylor series expansions in powers of IMF b z = B z /[ B 2

Fig. 10 .
Fig. 10.Equatorial plots of the electric current density, derived from the TS07D modelling of the storm of 21-23 April: early main phase (A), peak of the main phase (B), early recovery phase (C).

Fig. 11 .
Fig. 11.Comparing the actual variation of the SYM-H index (heavy black line) during the storm of 3-14 September, 2002, with that calculated from the TS05 model (thin black line).Separate contributions from individual current systems are shown by coloured lines in the upper panel.Centre and bottom panels show the concurrent variation of IMF Bz and the solar wind speed, respectively(Tsyganenko and Sitnov, 2005).

29 Fig. 11 .
Fig. 11.Comparing the actual variation of the SYM-H index (heavy black line) during the storm of 3-14 September, 2002, with that calculated from the TS05 model (thin black line).Separate contributions from individual current systems are shown by coloured lines in the upper panel.Centre and bottom panels show the concurrent variation of IMF B z and the solar wind speed, respectively(Tsyganenko and Sitnov, 2005).

Fig. 12 .
Fig. 12. Spatial coverage of the magnetosphere by space magnetometer data of the Polar (red), Geotail (light blue), Cluster (green), Themis-A, -D, -E (dark blue), and Themis-B, -C (yellow) spacecraft.The data set comprises 241,418 data records corresponding to 123 storm events between 1996 and 2012, with peak negative SYM-H index values not exceeding -200 nT.Equatorial and noon-midnight projections are shown in the left and right panels, respectively.

31Fig. 12 .
Fig. 12. Spatial coverage of the magnetosphere by space magnetometer data of the Polar (red), Geotail (light blue), Cluster (green), Themis-A, -D, -E (dark blue), and Themis-B, -C (yellow) spacecraft.The data set comprises 241 418 data records corresponding to 123 storm events between 1996 and 2012, with peak negative SYM-H index values not exceeding −200 nT.Equatorial and noon-midnight projections are shown in the left and right panels, respectively.near-perigee data are usually contaminated by large errors, rapidly increasing with decreasing r due to the fast-growing main geomagnetic field and inaccuracies of the satellite attitude data.The final step is to average the external field components over 5-minute intervals.To illustrate the overall coverage of the magnetosphere by space magnetometer data, Fig.12 displays the spatial distribution of observations in a grand modelling set containing 241,138 data records, used in the calibration of the most recent T13 model.The set includes 123 storm events during the period 1997-2012, with peak negative values of SYM-H not exceeding -200 nT.

Fig. 13 .
Fig. 13.Radial distributions of the data point density in the modelling set shown in Fig.12 (top panel), average B magnitude (centre panel), and the weight function, introduced in the fitting of the T13 model to data (bottom).

Fig. 13 .
Fig. 13.Radial distributions of the data point density in the modelling set shown in Fig. 12 (top panel), average B magnitude (centre panel), and the weight function, introduced in the fitting of the T13 model to data (bottom).

Fig. 14 .
Fig. 14.Scatter plots of the model field components vs the observed components, based on the data set shown in Fig.12.Only data taken outside r = 6.6 RE are included.Plots in the top and bottom rows compare the components of the total and external parts of the field, respectively.The corresponding correlation coefficients and regression slopes are shown in each panel.

925Fig. 14 .
Fig. 14.Scatter plots of the model field components vs. the observed components, based on the data set shown in Fig. 12.Only data taken outside r = 6.6 R E are included.Plots in the top and bottom rows compare the components of the total and external parts of the field, respectively.The corresponding correlation coefficients and regression slopes are shown in each panel.
(s) − b m (s)] ds ,(36)    where b a = B a /B a and b m = B m /B m are the actual and the model field direction vectors, respectively.Since the actual field vectors, and hence their deviations from the model δb are known only at discrete locations, irregularly scattered over the modelling region, the integrand δb(s) in Eq. (36) was approximated by a simple analytical function of the distance s along the model field line L m , whose parameters were external field model (black).As can be seen from the plots, the most recent model yields the best result, though with only a token improvement in the distant magnetosphere.

Fig. 15 . 35 Fig. 15 .
Fig. 15.Histograms of the angular difference between the observed and model vectors of the total magnetic field, based on the multi-spacecraft modelling data set displayed in Fig.12.Red, blue, and black lines correspond to the T96, T13, and zero external model field.The left and right panels show results for the inner (R < 6.6 RE) and outer (R > 6.6 RE) magnetosphere, respectively.The corresponding median values are shown by the vertical dashed lines.
between the ionosphere (s = 0) and the equatorial plane (s = Seq), so that s) − bm(s)] ds , (36) where ba = Ba/Ba and bm = Bm/Bm are the actual and the model field direction vectors, respectively.Since the actual field vectors, and hence their deviations from the model δb are known only at discrete locations, irregularly scattered over the modelling region, the integrand δb(s) in Eq.(36) 955 was approximated by a simple analytical function of the distance s along the model field line Lm, whose parameters were defined using data points inside a field line tube of a finite thickness, centered on the field line Lm.The above method was used to test the T89 model, with results presented in the form of polar diagrams displaying 2D distributions of the mapping errors for several levels of the Kp-index (see Plates 1-4 in the above-cited paper).The field line mapping direction can in 960 principle be reversed, so that one could evaluate the field line footpoint deviations in the ionosphere, instead of at their apex points.

Fig. 16 .Fig. 16 .
Fig. 16.Illustrating the IB concept.Left: a 3D view of a numerically traced particle orbit, crossing the equatorial current sheet outside of IB.After several bounces, the particle ends up in the ionosphere due to the non-adiabatic scattering of its magnetic moment into the loss cone.Right: trapped (red) and precipitating (blue) particle flux variation as seen by a low-altitude polar-orbiting satellite; poleward from the IB (green) both fluxes become equal.Another powerful technique to assess the models' performance in terms of the mapping accuracy is based on low-altitude observations of energetic charged particles.The idea of the method reduces to the simple fact that, for each species of particles with a given rigidity mV c/q there exists a surface 965

Fig. 17 .
Fig. 17.Radial profiles of the equatorial perpendicular (black solid line) and parallel (dash-dotted line) pressures, obtained by minimizing the rms difference between the magnetic and plasma stresses in the T96 model (a quiet-time case).For comparison, isotropic pressure profiles by Lui and Hamilton (1992) (blue), Tsyganenko and Mukai (2003) (green), and that based on 1D integration of the same T96 field (red) are shown (from Tsyganenko, 2010).

Fig. 17 .
Fig. 17.Radial profiles of the equatorial perpendicular (black solid line) and parallel (dash-dotted line) pressures, obtained by minimizing the rms difference between the magnetic and plasma stresses in the T96 model (a quiet-time case).For comparison, isotropic pressure profiles by Lui and Hamilton (1992) (blue), Tsyganenko and Mukai (2003) (green), and that based on 1-D integration of the same T96 field (red) are shown (from Tsyganenko, 2010).

Fig. 18 .
Fig. 18.Illustrating the effect of the magnetopause flaring rate on the distant tail configuration in the T13 model.Larger flaring due to negative IMF Bz results in more stretched field lines.

42 Fig. 18 .
Fig. 18.Illustrating the effect of the magnetopause flaring rate on the distant tail configuration in the T13 model.Larger flaring due to negative IMF B z results in more stretched field lines.

Fig. 19 .
Fig. 19.The top panel shows the variation of the total R1 (dark purple) and R2 (green) FACs (in MA) during the moderate storm of 8-10 March, 2008, as reproduced by the T13 model.The centre and bottom panels show the concurrent variations of IMF B z and the SYM-H index, respectively.

Fig. 20 .Fig. 20 .
Fig. 20.Polar cusps in an idealized vacuum model with only a pair of null B points (left), and in the real magnetosphere, with deep intrusions of the magnetosheath plasma (right).1165The second important phenomenon associated with the polar cusps is a high correlation of the B ycomponent inside the cusp funnels with IMF B y .As found in a statistical study of that effect based on Polar and Cluster data(Tsyganenko, 2009), the linear regression coefficient relating B cusp y with