The adjoint sensitivity method of global electromagnetic induction for CHAMP magnetic data

An existing time-domain spectral-ﬁnite element approach for the forward modelling of electromagnetic induction vector data as measured by the CHAMP satellite is, in this paper, supplemented by a new method of computing the sensitivity of the CHAMP electromagnetic induction data to the Earth’s mantle electrical conductivity, which we term the adjoint sensitivity method. The forward and adjoint initial boundary-value problems, both solved in the time domain, are identical, except for the speciﬁcation of prescribed boundary conditions. The respective boundary-value data at the satellite’s altitude are the X magnetic component measured by the CHAMP vector magnetometer along the satellite track for the forward method and the difference between the measured and predicted Z magnetic component for the adjoint method.The squares of these differences summed up over all CHAMP tracks determine the misﬁt. The sensitivities of the CHAMP data, that is the partial derivatives of the misﬁt with respect to mantle conductivity parameters, are then obtained by the scalar product of the forward and adjoint solutions, multiplied by the gradient of the conductivity and integrated over all CHAMP tracks. Such exactly determined sensitivities are checked against numerical differentiation of the misﬁt, and good agreement is obtained. The attractiveness of the adjoint method lies in the fact that the adjoint sensitivities are calculated for the price of only an additional forward calculation, regardless of the number of conductivity parameters. However, since the adjoint solution proceeds backwards in time, the forward solution must be stored at each time step, leading to memory requirements that are linear with respect to the number of steps undertaken. Having determined the sensitivities, we apply the conjugate gradient method to infer 1-D and 2-D conductivity structures of the Earth based on the CHAMP residual time series (after the subtraction of static ﬁeld and secular variations as described by the CHAOS model) for the year 2001. We show that this time series is capable of resolving both 1-D and 2-D structures in the upper mantle and the upper part of the lower mantle, while it is not sufﬁciently long to reliably resolve the conductivity structure in the lower part of the lower mantle. of to the allowing study of how latitudinal in to assimilation of at into


1374
Z. Martinec and J. Velímský (2) The satellite is assumed to orbit the Earth sufficiently fast (one CHAMP orbit takes approximately 90 min) compared to the time variations of the ring current. This assumption allows the separation of the spatial and temporal changes of the magnetic field observed by the single satellite to be performed in a simple way. Each night-side satellite track is considered to sample a snapshot of the magnetic field at the time when the satellite crosses the magnetic equator.
(3) Since the CHAMP satellite orbit is nearly polar (the inclination of the CHAMP orbit is approximately 87.3 • ), magnetic signals sampled along a track are dominantly influenced by latitudinal changes in the electrical conductivity of the Earth's mantle. We will not consider the effect of longitudinal variations in electrical conductivity on track data and assume that the electrical conductivity σ of the Earth is axially symmetric, that is where B is a conducting sphere approximating a heterogeneous Earth, r is the radial distance from the centre of B and ϑ is the colatitude.
(4) Since ring-current magnetospheric excitation has nearly an axially symmetric geometry, we will assume that, for a given satellite track, the inducing and induced magnetic fields possess axial symmetry, that is where G (e) jm (t) and G (i) jm (t) are the time-dependent, spherical harmonic Gauss expansion coefficients of the magnetic field generated by external equatorial ring currents and the magnetic field generated by the induced eddy currents in the Earth, respectively.

Forward method
The conducting sphere B is assumed to be surrounded by an insulating atmosphere, approximated by a spherical layer A with the inner boundary coinciding with the surface ∂B of the sphere B with the mean Earth radius r = a and the outer boundary with the mean-orbit sphere ∂A of radius r = b. The solution domain G for EM induction modelling is the unification of the conducting sphere B and the insulating spherical layer A, that is G = B ∪ A, bounded by the mean-orbit sphere ∂G = ∂A.
The magnetic field in G is induced by time-varying electrical currents in the magnetosphere with characteristic timescales ranging from several hours to tens of days. For the axisymmetric geometry of external ring-current sources and conductivity models, it is convenient to formulate the initial, boundary-value problem (IBVP) of global EM induction in terms of the toroidal vector potential A such that the magnetic induction vector is expressed in the form B = curl A. The mathematical formulation is as follows: find the toroidal vector potential A such that with the boundary condition where the magnetic permeability μ of G is assumed to be constant. In this paper, eqs (3) and (4) are applied for times t ∈ (0, T ). The term B t represents the tangential component of the magnetic induction vector B at the satellite altitude and n is the unit normal to ∂G. The axisymmetric geometry allows us (see Section 4) to determine B t from the horizontal northward X component of the magnetic induction vector B measured by the CHAMP magnetometer. Moreover, the toroidal vector potential A is considered divergence-free, that is div A = 0. Note that the conductivity σ = 0 in the insulating atmosphere A implies that the second term in eq. (3) vanishes in A. Eq. (3) is yet subject to the initial condition where A 0 is the generating potential for the initial magnetic induction B 0 such that B 0 = curl A 0 .

A D J O I N T S E N S I T I V I T Y M E T H O D O F E M I N D U C T I O N F O R T H E Z C O M P O N E N T O F C H A M P M A G N E T I C DATA
In this section, we formulate the adjoint sensitivity method of EM induction for computing the sensitivity of the Z component of CHAMP magnetic data with respect to mantle-conductivity structure.

Misfit function and its gradient in the parameter space
Let us first consider the conductivity σ (r , ϑ) of the conducting sphere B to be represented in terms of an M-dimensional system of rand ϑ-dependent base functions and let the expansion coefficients of this representation be σ 1 , σ 2 , . . . , σ M . Defining the parameter vector σ := (σ 1 , σ 2 , . . . , σ M ), the dependence of the conductivity σ (r , ϑ) on the conductivity parameters σ can be made explicit as σ = σ (r, ϑ; σ ).
The adjoint method for EV induction 1375 Martinec & McCreadie (2004) showed that the solution of the IBVP for CHAMP magnetic data formulated in the previous section enables the modelling of the time evolution of the normal component B n := n · B of the magnetic induction vector on the mean-orbit sphere ∂ G along the satellite tracks. These predicted data B n ( σ ) can be compared with the observations B (obs) n of the normal component of the magnetic induction vector by the CHAMP onboard magnetometer. The differences between observed and predicted values can be used as a misfit for the inverse EM-induction modelling. The adjoint method of EM induction presented in this paper calculates the sensitivity of the forward-modelled data B n ( σ ) on the conductivity distribution σ by making use of the differences B (obs) n − B n ( σ ) as boundary-value data.
Let the observations B (obs) n be made for times t ∈ (0, T ) such that, according to assumption (ii), B (obs) n (ϑ, t i ) at a particular time instance t i ∈ (0, T ) corresponds to the CHAMP observations along the ith satellite track. The least-squares misfit is therefore defined as where the weighting factor w b = w b (ϑ, t) is chosen to be dimensionless such that the misfit has the SI unit m 3 sT 2 /[μ], [μ] = kg m s −2 A −2 . If the observations B (obs) n contain random errors, which are statistically independent, the statistical variance of observations may be substituted for the reciprocal value of w 2 b (e.g. Bevington 1969, section 6-4). In Section 4.2, the ϑ dependence of w b allows us to eliminate the track data from the polar regions which are contaminated by signals from field-aligned currents and polar electrojets while the time dependence of w b allows us to eliminate the track data for time instances when other undesirable magnetic effects at low and mid-latitudes contaminate the signal excited by equatorial ring currents.
The sensitivity analysis or inverse modelling requires the computation of the partial derivative of the misfit with respect to the model parameters, that is the derivatives ∂χ 2 /∂ σ m , m = 1, . . . , M, often termed the sensitivities of the misfit with respect to the model parameters σ m (e.g. Sandu et al. 2003). To abbreviate the notation, the partial derivatives with respect to the conductivity parameters are ordered in the gradient operator in the M-dimensional parameter space, where the hat inσ m indicates a unit vector.
Realizing that the observations B (obs) n are independent of the conductivity parameters σ , that is ∇ σ B (obs) where B n ( σ ) are the weighted residuals of the normal component of magnetic induction vector, that is, The straightforward approach to find ∇ σ χ 2 , is to approximate ∂χ 2 /∂σ m by a numerical differentiation of forward model runs. Due to the size of the parameter space, this procedure is often extremely computationally expensive.

The forward sensitivity equations
The forward sensitivity analysis computes the sensitivities of the forward solution with respect to the conductivity parameters, that is the partial derivatives ∂ A/∂σ m , m = 1, . . . , M. Using them, the forward sensitivities ∇ σ B n are computed and substituted into eq. (9) for ∇ σ χ 2 .
To form the forward sensitivity equations, also called the linear tangent equations of the model (e.g. McGillivray et al. 1994;Cacuci 2003;Sandu et al. 2003Sandu et al. , 2005, let us consider the conductivity model (6) in eq. (3) and differentiate the forward model eqs (3)-(5) with respect to the conductivity parameters σ with homogeneous boundary and initial conditions and where ∇ σ B t = ∇ σ A 0 = 0 have been substituted because the boundary data B t and the initial condition A 0 are independent of the conductivity parameters σ . In the forward sensitivity analysis, for each parameter σ m and associated forward solution A, a new source term ∇ σ σ ∂ A/∂t is created and the forward sensitivity eqs (11)-(13) are solved to compute the partial derivative ∂ A/∂σ m . The forward sensitivity analysis is known to be very effective when the sensitivities of a larger number of output variables are computed with respect to a small number of input parameters (Sandu et al. 2003;Petzold et al. 2006).

The adjoint sensitivity equations
The adjoint method provides an efficient alternative to the forward sensitivity analysis for evaluating ∇ σ χ 2 without explicit knowledge of ∇ σ A, that is without solving the forward sensitivity equations. Hence, the adjoint method is more efficient for problems involving a large number of model parameters.
The adjoint sensitivity analysis proceeds by forming the inner product of eqs (11) and (12) with a yet unspecified adjoint functionÂ. Because the forward sensitivity equations are linear in ∇ σ A, an adjoint equation exists (Cacuci 2003). To derive it, eqs (11) and (12) are multiplied by an adjoint functionÂ(r, ϑ, t), then integrated over G and ∂ G, respectively, and subtracted from each other: where dot stands for the scalar product of vectors.
In the next step, the integrals in eq. (14) are transformed such that ∇ σ A andÂ interchange. To achieve this, let us consider the Green's theorem for two sufficiently smooth functions f and g in the form Interchanging the functions f and g and subtracting the new equation from the original one results in the integral identity By this, the positions of ∇ σ A andÂ can be exchanged in the first two integrals in eq. (14): To perform the same transformation in the third integral, we integrate eq. (17) over the time interval t ∈ (0, T ), that is Then we exchange the order of integration over the spatial variables and time in the third integral and perform the time integration by partes: The second term on the right-hand side is equal to zero because of homogeneous initial condition given by eq. (13). Finally, eq. (18) takes the form Remembering that ∇ σ B n is the derivative that we wish to eliminate from ∇ σ χ 2 , we add the homogeneous eq. (20) to eq. (9) (note the physical units of eq. (20) are the same as ∇ σ χ 2 , namely m 3 sT 2 /[μσ ] provided that the physical units ofÂ are the same as of A, namely, Tm): The adjoint functionÂ has been considered arbitrary so far. Our aim is now to impose constraints on it such that the originally arbitrarŷ A transforms to the well-defined adjoint toroidal vector potential. We first eliminate the volume integrals over G proportional to ∇ σ A by requiring that with the terminal condition onÂ: The boundary condition forÂ on ∂G is derived from the requirement that the surface integrals over ∂G in eq. (21) cancel each other, that is The adjoint method for EV induction 1377 at any time t ∈ (0, T ). We will elaborate on this condition in the next section. Under these constraints, the gradient of χ 2 ( σ ) takes the form

Boundary condition for the adjoint potential
To relate ∇ σ A and ∇ σ B n in the constraint described by eq. (24) and, subsequently, to eliminate ∇ σ A from it, a parameterization of A is necessary. For global EM induction, it is convenient to parameterize A in terms of vector spherical harmonics. Such a parameterization is described in detail by Martinec (1997) and Martinec et al. (2003). Here, we only introduce the final form of the representation of A, which is also employed for the adjoint potentialÂ: where Y j j (ϑ) are the zonal toroidal vector spherical harmonics. Their definition and some properties are given in the Appendix A. We first aim to express the gradient ∇ σ χ 2 in terms of spherical harmonics. Since the upper boundary ∂G of the solution domain G is the mean-orbit sphere of radius b, the external normal n to ∂G coincides with the unit vector e r , that is n = e r . Applying the gradient operator ∇ σ on the equation B n = e r · curl A and using eq. (A14), we obtain Moreover, applying a two-step least-squares analysis (Martinec & McCreadie 2004) on the residual satellite-track data B n defined by eq. (10), these observables can, at a particular time t ∈ (0, T ), be represented as a series of the zonal scalar spherical harmonics with spherical harmonic coefficients of the form Substituting eqs (27) and (28) into eq. (9) and employing the orthonormality property (A2) of the zonal scalar spherical harmonics Y j (ϑ), the gradient of the misfit χ 2 becomes We now return to the constraint (24) and express it in terms of spherical harmonics. By the parameterization (26) and the assumption n = e r , the differential relation (A16) applied toÂ yields where The first constituent in the first integral of the constraint (24) is expressed by eq. (31), while the second constituent can be obtained by applying the gradient operator ∇ σ to eq. (26). The two constituents in the second integral of the constraint (24) are expressed by eqs (27) and (28), respectively. Performing all indicated substitutions, we obtain Interchanging the order of integration over the full solid angle and summations over j's, and making use of the orthonormality properties (A2) and (A7) of the zonal scalar and vector spherical harmonics, respectively, eq. (33) reduces to Z. Martinec and J. Velímský which is to be valid at any time t ∈ (0, T ). To satisfy this constraint independently of ∇ σ A j j (b, t), one last condition is imposed upon the adjoint potentialÂ, namely at any time t ∈ (0, T ).

Adjoint method
We can now summarize the formulation of the adjoint method of EM induction for CHAMP satellite magnetic data.
Given the electrical conductivity model σ (r , ϑ) in the sphere B, the forward solution A(r , ϑ, t) in B and the atmosphere A for t ∈ (0, T ) and the observations B (obs) n (t) on the mean-orbit sphere ∂ G of radius r = b with uncertainties quantified by weighting factor w b , find the adjoint potentialÂ(r, ϑ, t) in G = B ∪ A by solving the adjoint problem: with the boundary condition and the terminal condition The gradient of the misfit χ 2 ( σ ) is then expressed as The set of eqs (36)-(38) is referred to as the adjoint problem of the forward problem specified by eqs (3)-(5). Combining the forward solution A and the adjoint solutionÂ according to eq. (39) thus gives the exact derivative of the misfit χ 2 .

Reverse time
The numerical solution of eq. (36), solved backwards in time from t = T to t = 0, is inherently unstable. Unlike the case of the forward model equation and the forward sensitivity equation, the adjoint equation effectively includes negative diffusion, which enhances numerical perturbations instead of smoothing them, leading to an unstable solution. To avoid such numerical instability, we change the sign of the diffusive term in eq. (36) by reversing the time variable. We introduce the reverse time τ = T − t, τ ∈ (0, T ), and the reverse-time adjoint potentialǍ(τ ) such that Hence, and eq. (36) transforms to the diffusion equation for the reverse-time adjoint potentialǍ(τ ): with the boundary condition The terminal condition (38) forÂ becomes the initial condition for the potentialǍ: With these changes, the adjoint equations become similar to those of the forward method, and thus nearly identical numerical methods can be applied. In Appendix B, we demonstrate the way of how the adjoint method forǍ can be reformulated in a weak sense. The gradient (39) transforms to The importance of eq. (45) is that, once the forward problem (3)-(5) is solved and the misfit χ 2 is evaluated from eq. (7), the gradient ∇ σ χ 2 may be evaluated for little more than the cost of a single solution of the adjoint system (42)-(44) and a single scalar product in eq. (45), regardless of the dimension of the conductivity vector σ . This is compared to other methods of evaluating ∇ σ χ 2 that typically require the solution of the forward problem (3)-(5) per component of σ . In Appendix C, the method for numerical evaluation of misfit gradient ∇ σ χ 2 is shown. We now explain the specific steps involved in the adjoint computations. First, the forward solutions A (t i ) are calculated at discrete time levels 0 = t 0 < t 1 < · · · < t n = T by solving the forward problem (3)-(5), and each solution A (t i ) must be stored. Then the reverse-time adjoint solutionsǍ(t i ), i = 0, . . . , n, are calculated, proceeding again forwards in time according to eqs (42)-(44). As each adjoint solution is computed, the misfit and its derivative are updated according to eq. (7) and (45), respectively. WhenǍ(T ) has finally been calculated, both χ 2 and ∇ σ χ 2 are known. The forward solutions A (t i ) are stored because eqs (43) and (45) depend on them for the adjoint calculation. As a result, the numerical algorithm has memory requirements that are linear to the number of timesteps. This is the main drawback of the adjoint method.

Spherical harmonic parameterization of CHAMP magnetic data
The spherical harmonic representation of the tangential component B t of magnetic data at satellite altitudes, used in the forward method of EM induction, can be obtained by substituting eq. (B2) into eq. (A16): where are the spherical harmonic coefficients of the horizontal northward X component of the magnetic induction vector B measured at satellite Similarly, the spherical harmonic representation of the normal component B n of magnetic data at satellite altitudes, used in the adjoint method of EM induction, can be obtained by substituting eq. (B2) into eq. (A14): where are the spherical harmonic coefficients of the Z component of the magnetic induction vector B measured at satellite altitude r = b, that is Eqs (47) and (50) show that the coefficients X j (t) and Z j (t) are composed of two different linear combinations of the spherical harmonics G (e) j (t) of the external EM sources and the spherical harmonics G (i) j (t) of the induced magnetic field inside the Earth. Consequently, there is no need to specify these coefficients separately when X j (t) and Z j (t) are used as the boundary-value data for the forward and adjoint modelling of EM induction, respectively.

Selection and processing of vector data
The data analyzed in this paper were recorded by a three-component vector magnetometer on board of CHAMP. To demonstrate the performance of the adjoint sensitivity method, we have selected, from all records spanning more than 8 years, the 1-year long time series from January 1, 2001 (track no. 2610), to January 10, 2002 (track no. 8402). Judging from the Dst index (Fig. 2), there were about 10 events when the geomagnetic field was significantly disturbed by geomagnetic storms or substorms. In order to minimize the effect of strong day-side ionospheric currents, we use only night-side data recorded by the satellite between 18:00 and 6:00 local solar time.
In the first step of the data processing, we use the CHAOS model of the Earth's magnetic field  to isolate the signals corresponding to EM induction by storm-time magnetospheric currents. Based on the CHAOS model, we remove the main and crustal fields up to degree 50 and the secular variation up to degree 18 from the CHAMP data. Left panels: the original CHAMP data plotted along geographical colatitude. X g , Y g and Z components point, respectively, to the geographic north, the geographic east and downwards. Right panels: Black lines denote X -and Z-CHAMP components after the removal of the CHAOS model and the rotation of the residual field to dipole coordinates. The red lines show the results of the two-step, track-by-track spherical harmonic analysis, including the extrapolation into the polar regions by using data from the mid-colatitude interval (40 • , 140 • ), as marked by dotted lines.
In the next step, the horizontal magnetic components (X , Y ) are rotated from geographic coordinates to dipole coordinates, assuming that the north geomagnetic pole is at 78.8 • N, 70.7 • W. Since we assume an axisymmetric geometry of external currents and mantle electrical conductivity, the dipolar longitudinal component Y is not considered hereafter and we use X and Z to describe, respectively, the northward and downward magnetic components in dipolar coordinates. Fig. 1 shows an example of the original and processed data from CHAMP track no. 6753.

Spherical harmonic analysis
The two-step, track-by-track spherical harmonic analysis (Martinec & McCreadie 2004) is applied to the X component of the CHAMP track data, resulting in the least-squares estimate X (obs) j (t) of the coefficients X j (t). This method allows us to exclude observations from the polar regions, which are contaminated by signals from field-aligned currents and polar electrojets. Instead, the field in these regions is extrapolated from the field at low and mid-latitudes in accordance with the assumption that global EM induction is driven by the equatorial ring current. The extrapolation takes advantage of the fact that the parameterization (48) ensures the X component approaches zero at the northern and southern magnetic poles.
The crucial points of the extrapolation are the selection of the truncation degree j max of the parameterization (48) and the determination of the colatitude interval (ϑ 1 , ϑ 2 ) where the data are not disturbed by the polar currents. Velímský et al. (2006) imposed three criteria to determine these two parameters. First, the power of the magnetic field from the external ring currents is concentrated in low-degree harmonic coefficients, particularly in the j = 1 term, and the leakage of EM energy into higher-degree terms caused by the Earth's conductivity and electrical-current geometry is monotonically decreasing. Second, the first derivative of the X component with respect to colatitude does not change sign in the polar regions to exclude unrealistic oscillatory behavior of the X component in these regions caused by a high-degree extrapolation. Third, if the least-squares estimate of the X component of CHAMP data over the colatitude interval (ϑ 1 − 5 • , ϑ 2 + 5 • ) differs by more than 10 nT compared to the estimate over the interval (ϑ 1 , ϑ 2 ), the field due to the polar currents is assumed to encroach upon the field produced by near-equatorial currents, and the narrower colatitude interval (ϑ 1 , ϑ 2 ) is considered to contain only the signature generated by the near-equatorial currents. Applying these criteria to the CHAMP track data iteratively, starting from degree j = 1 and the colatitude interval (10 • , 170 • ) and proceeding to higher degrees and shorter colatitude intervals, we found that the maximum cut-off degree is equal to j max = 4 and the colatitude interval is (40 • , 140 • ).
The same two-step, track-by-track analysis is applied to the Z component of CHAMP track data, providing the least-squares estimates Z (obs) j (t) of the coefficients Z j (t). However, the extrapolation of the Z component from the field at low and mid-latitudes is more problematic (t) (blue) of horizontal and vertical components obtained by the two-step, track-by-track spherical harmonic analysis of CHAMP data for the year 2001. A mean and linear trend are removed following the arguments of Olsen et al.(2005). The coefficients from the missing tracks are filled by cubic spline interpolation applied to the detrended time series. Note that the sign of the X 1 component is opposite to that of the Dst index (black line). This is because this component is expressed in terms of the associated Legendre functions with the norm by Varshalovich et al. (1989), which differs from the Schmidt quasi-normalization, except a different normalization, by the factor (−1) m , where m is the azimuthal order of the associated Legendre functions. Time on the horizontal axis is measured from midnight of January 1, 2000. than that for the X component. This is because (1) the second selection criterion cannot be applied since the Z component does not approach zero at the magnetic poles as seen from parameterization eq. (51), and (2) the Z component of CHAMP magnetic data contains a larger portion of high-frequency noise than the X component, which, in principle, violates the assumption of the third selection criterion. Fig. 3 shows that the leakage of EM energy from the j = 1 into higher-degree terms is not monotonically decreasing for the Z component. That is why the least-squares estimates Z (obs) j (t) are extrapolated to polar regions from the colatitude interval (ϑ 1 , ϑ 2 ) and up to the spherical degree j max determined for the X component.
The procedure applied to the 2001-CHAMP track data results in the time series of spherical harmonic coefficients X (obs) j (t) and Z (obs) j (t), j = 1, . . . , 4. As an example, the resulting coefficients for degree j = 1 are plotted in Fig. 2 as functions of time after January 1, 2001. As expected, there is a high correlation between the first-degree harmonics X  (t) (bottom panel) components. Degrees of 1, 2, 3 and 4 are, respectively, shown by black, red, blue and green lines. The spectra have peaks at higher harmonics of the 27-day solar rotation period, that is at periods of 9 days, 6.8 days, 5.6 days, 4.8 days etc.

Power-spectrum analysis
Although the method applied in this paper is based on the time-domain approach, it is valuable to inspect the spectra of the X (obs) j (t) and Z (obs) j (t) time series. Fig. 3 shows the maximum-entropy power-spectrum estimates (Press et al. 1992, section 13.7) of the first four spherical harmonics of the horizontal and vertical components. We can see that the magnitudes of the power spectra of the X component monotonically decrease with increasing harmonic degree, which is a consequence of the first selection criterion applied in the two-step, track-by-track analysis. For instance, the power spectrum of the second-degree terms are about two orders of magnitude smaller than that of the first-degree term. As already introduced, and also seen in the bottom panel of Fig. 3, this is not the case for the Z component. The magnitude of the maximum-entropy power-spectrum of the Z component is larger than that of the X component for j > 1, which demonstrates that the Z component of CHAMP magnetic data contains a larger portion of high-frequency noise than the X component.
Despite analyzing only night-side tracks, there is a significant peak at the period of 1 day in the power spectra of higher degree harmonics ( j ≥ 2), but, surprisingly, missing in the spectra of the first-degree harmonic. To eliminate the induction effect of residual dawn/dusk ionospheric electrical currents, we shrink the night-side local-solar time interval from (18:00, 6:00) to (22:00, 4:00). However, a 1-day period signal remains present in the CHAMP residual signal (not shown here). To locate a region of potential inducing electrical currents, time series of X (obs) j (t) and Z (obs) j (t) coefficients are converted to time series of spherical harmonic coefficients of external and internal fields counted with regard to the CHAMP satellite altitude. To obtain these coefficients, denoted byĜ (e) j andĜ (i) j , the Gaussian expansion of the external magnetic potential is taken at satellite orbit of radius r = b, which results in eqs (47) and (50) where the mean Earth's radius a is replaced by b. The straightforward derivation then yieldŝ The maximum-entropy power-spectrum estimates of external and internal coefficientsĜ (e) j andĜ (i) j are shown in Fig. 4. We can see that these spectra for degrees j = 2 to 4 have also a peak at a period of 1 day. This means that at least part ofĜ (e) j originates in the magnetosphere or  even magnetopause and magnetic tail, whileĜ (i) j may originate from the residual night-side ionospheric currents and/or electrical currents in the Earth induced by either effect. The 1-day period signal may also result from the fact that we express the magnetospheric field in the dipolar coordinate system while 1-day periodicities may be better represented in the solar magnetospheric coordinates related to the position of the Sun (Kivelson & Rusell 1995;Maus & Lühr 2005).
In addition, Fig. 4 shows that, while the periods of peak values in the external and internal magnetic fields for degree j = 1 correspond to each other, for the higher degree harmonic terms, such a correspondence is only valid for some periods, for instance 6.8, 5.6 or 4.8 days. However, the peak for the period of 8.5 days in the internal component for j = 2 is hardly detectable in the external field. This could be explained by a three-dimensionality effect in the electrical conductivity of the Earth's mantle that causes the leakage of EM energy from degree j = 1 to the second and higher degree terms. This leakage may partly shift the characteristic periods in the resulting signal due to interference between signals with various spatial wavelengths and periods.

V E R I F I C AT I O N
Before the adjoint sensitivities are computed for the CHAMP data, we will, as a first step, verify the adjoint sensitivity method (ASM) against the direct numerical differentiation of eq. (7) for a simple three-layer, 1-D conductivity model.

Brute-force sensitivities
Sensitivities generated with the adjoint sensitivity method, called hereafter as the adjoint sensitivities, are now compared to those generated by direct numerical differentiation of the misfit, the so-called brute-force method (BFM) (e.g. Bevington 1969), in which the partial derivative of misfit with respect to σ m at the point σ 0 is approximated by the second-order accuracy, centred difference of two forward model runs: where ε refers to a perturbation applied to the nominal value of σ 0 m .

Model parameterization
According to eq. (C1), the electrical conductivity is parameterized by the functions σ (ϑ), which vary only with co-latitude ϑ within the interval R ≤ r ≤ R +1 , = 1, · · ·, L. Moreover, we parameterize σ (ϑ) by the zonal scalar spherical harmonics Y j (ϑ). As a result, the logarithm of the electrical conductivity is considered in the form where ξ (r ) is equal to 1 in the interval R ≤ r ≤ R +1 ; 0 elsewhere. The number of conductivity parameters σ j , that is the size of conductivity parameter vector σ , is M = L(J + 1). Martinec et al. (2003) developed a time-domain spectral-finite element method (the TISFEM method) for ground magnetic data. This method is governed by the variational equation

Three-layer, 1-D conductivity model
where the sesquilinear form a 1 (·, ·) and the functional F 1 (·, ·) are specified in Martinec et al. (2003). They tested the TISFEM method by comparing it with the analytical and semi-analytical solutions of EM induction for the cases of two concentrically and eccentrically nested spheres of different, but constant, electrical conductivities. They showed that the numerical code implementing the TISFEM method for ground magnetic data performs correctly, and that the TISFEM method is particularly appropriate when the external current excitation is transient. Martinec & McCreadie (2004) modified the TISFEM method for satellite magnetic data and used it for ground magnetic data to verify this modification. We now use both the original and modified TISFEM methods to verify the ASM for computing the partial derivatives of the misfit with respect to the conductivity parameters σ 0 . The complex structure of a magnetic storm will be described by a simple mathematical model that simulates a storm's basic features. The storm ring current is considered axisymmetric with a P 10 (cos ϑ) spatial structure. Consequently, all spherical harmonics of the external scalar magnetic potential are equal to zero, except for the first-degree coefficient G (e) 1 (t). After the onset of a magnetic storm at t = 0, the ring current quickly peaks and then decays exponentially. This time evolution is modelled by the function (Martinec et al. 2003) where √ 4π/3 is the inverse norm of P 10 (cos ϑ), A is the amplitude and τ is the relaxation time describing the recovery phase of the storm. We will use τ = 1 day and A = 0.003 nT s −1 in the following test example.
Consider a 1-D conducting sphere B, consisting of the lithosphere, the upper mantle (UM), the upper (ULM) and lower (LLM) parts of the lower mantle, and the core. The interfaces between the conductivity layers are kept fixed at depths of 220, 670, 1500 and 2890 km, respectively. The conductivities of the lithosphere and the core are 0.001 S/m and 10000 S/m, respectively, and fixed at these values for all sensitivity tests, hence the number of conductivity parameters σ 0 is L = 3. The nominal values of the conductivity parameters are σ 0 10 = 1 (hence σ LLM = 10 S/m), σ 0 20 = 0 (σ ULM = 1 S/m), and σ 0 30 = −1 (σ UM = 0.1 S/m). The strategy for verifying the ASM for satellite magnetic data consists of the following steps. First, the conducting sphere B is excited by a source field with the time evolution (56) applied as the boundary-value data on the surface ∂B. The numerical code implementing the TISFEM method for ground magnetic data, which is governed by the variational eq. (55), then computes the EM induction response of B. This provides the coefficients G (i) j (t), j = 1, 2, · · ·, of the induced magnetic field. Second, the external and internal coefficients are used to generate the boundary-value data B t at the satellite altitude according to eqs (46)-(47). These data are applied on the mean-orbit sphere ∂A and excite an electromagnetic field in the model consisting of the conducting sphere B and non-conducting atmosphere A. The numerical code implementing the forward TISFEM method for the X component of the satellite magnetic data, which is governed by the variational eq. (B7) with the boundary-value data (B15) 1 , computes the Z component of the EM induction response of the model G = B ∪ A. Thirdly, the conductivity parameters of the conducting sphere B are changed and the numerical code implementing the forward TISFEM method for the X component of satellite magnetic data is again applied, yielding the Z component of the EM induction response at satellite altitude. Since the conductivity structure of B has been changed, the new Z data differ from those computed in the previous step for the structure of B with the nominal conductivities. Fourth, the differences in the Z component, computed according to eq. (B15) 2 , are applied as the boundary-value data in the ASM and the adjoint sensitivities are computed. Finally, the derivatives of the misfit obtained by the ASM are compared with those computed by the BFM. Note that the sensitivity tests are carried out in such a way that the conductivity parameter of only one layer is varied at a time, with the conductivities of the other two layers kept fixed and equal to the nominal values.
Both methods used to compute the partial derivatives of the misfit with respect to the conductivity parameters suffer from approximation errors, which we will now compare. The approximation error of the BFM, which is based on the second-order differences (53), is proportional to the third-order derivative of the misfit multiplied by the square of the change ε in the conductivity parameter σ 0 . Consequently, the smaller the value of ε, the smaller the approximation error of the BFM.   Fig. 5, but for the approximation error of the adjoint sensitivities ∇ σ χ 2 for various time steps t (blue diamonds stand for t = 90 min, green plus signs for t = 45 min and red crosses for t = 22.5 min). The black dots show the brute-force sensitivities computed for a fixed value of ε = 0.01. Fig. 5 shows a comparison between the brute-force sensitivities for various values of ε and the adjoint sensitivities. The forward and adjoint solutions are computed for the mathematical model (56) of a storm with a time step t = 1.5 h and the three-layer, 1-D conductivity model described previously. We vary the size of ε in the conductivity parameters σ 0 used in the BFM while the time step t used in the ASM is kept fixed. Hence, the approximation error of the BFM is changing, while for the ASM it is fixed. Fig. 5 shows that the approximation error of the BFM decreases when ε decreases, however, there is a threshold of ε (≈0.01, black crosses in Fig. 5), below which, the differences between the BFM and ASM do not decrease in magnitude (for instance, for ε = 0.001, red dots in Fig. 5). This indicates that for this case, the approximation error of the ASM prevails. Martinec et al. (2003) approximated the time derivative of the toroidal vector potential, that is the term ∂ A/∂t which occurs in the variational equality (B7), and the gradient of the misfit (C2) by the forward Euler differencing scheme which introduces the approximation error proportional to the second-order derivatives of A multiplied by t. This error, which propagates into the forward and adjoint solutions, can be reduced by decreasing the size of either t (Fig. 6), or ∂ 2 A/∂t 2 (Fig. 13). Fig. 6 compares the adjoint sensitivities for various sizes of t with the brute-force sensitivities. The size of ε for the later ones is kept fixed, equal to the threshold value taken from the previous case study. Fig. 6 demonstrates that the approximation error of the ASM decreases with decreasing t, and can be reduced to a level such that the differences between the adjoint and brute-force sensitivities are about two orders of magnitude smaller than the sensitivities themselves. For the 2001-CHAMP magnetic data used in the next section, we will choose ε = 0.01 and t = 1 h.

S E N S I T I V I T Y A N A LY S I S F O R C H A M P DATA
The forward and adjoint solutions are now computed for the 2001-CHAMP data (see Section 4) with spherical-harmonic cut-off degree j max = 4 and time step t = 1 h. We will perform the sensitivity analysis of the data with respect to two different conductivity models, three-layer, 1-D conductivity model and two-layer, 2-D conductivity model. For each case, we first investigate the approximation error of the ASM and then run the conjugate gradient inversion for searching an optimal conductivity model adjusting the Z component of CHAMP data in least squares sense.

Sensitivity comparison
The first check of the adjoint sensitivities is performed for the three-layer, 1-D conductivity model described in Section 5.3. The results of the sensitivity tests are summarized in Fig. 7, where the top panels show the misfit χ 2 as a function of one conductivity parameter σ 0 whereas the other two are equal to the nominal values. The bottom panels compare the derivatives of the misfit obtained by the ASM with the BFM. From these results, two conclusions can be drawn. First, the differences between the derivatives of the misfit obtained by the ASM and BFM (the dashed lines in the bottom panels) are about one order (for σ 30 ) and at least two orders (for σ 10 and σ 20 ) of magnitude smaller than the derivatives themselves. Recalling the results of the verification test shown in Fig. 6, it justifies the validity of the ASM. The differences between the adjoint and brute-force sensitivities are caused by the approximation error of time numerical differentiation (57). In Section 6.3., we will show the way of how this error can further be reduced. Second, both the top and bottom panels show that the misfit χ 2 is most sensitive to the conductivity changes in the upper mantle and decreases with increasing depth of the conductivity layer, being least sensitive to conductivity changes in the lower part of the lower mantle.

Conjugate gradient inversion
The sensitivity results in Fig. 7 are encouraging with regard to the solution of the inverse problem for a 1-D mantle conductivity structure. We employ the conjugate gradient (CG) minimization with bracketing and line searching using Brent's method with derivatives (Press et al. 1992, section 10.3) obtained by the ASM. The inverse problem is solved for the 3 parameters σ 0 , with starting values equal to (1.5, 0, −1). Fig. 8 shows the results of the inversion, where the left panel displays the conductivity structure in the three-layer mantle and the right panel the misfit χ 2 as a function of the CG iterations. The blue line shows the starting model of the CG minimization, the dotted line the model after the first iteration and the red line the model after 10 iterations. As expected from the sensitivity tests, the minimization first modifies the conductivities of the UM and ULM, to which the misfit χ 2 is the most sensitive. When the UM and ULM conductivities are improved, the CG minimization also changes the LLM conductivity. The optimal values of the conductivity parameters after 10 iterations are (σ 10 , σ 20 , σ 30 ) = (1.990, 0.186, − 0.501). This corresponds to the conductivities σ ULM = 1.53 S/m and σ UM = 0.32 S/m for ULM and UM, which are considered to be well resolved, while the conductivity σ LLM = 97.8 S/m should be treated with some reservation, because of its poor resolution. A CHAMP time series longer than 1 year would be necessary to increase the sensitivity of CHAMP data to the LLM conductivity.

Sensitivity comparison
The second check of the adjoint sensitivities is performed for the 2-D conductivity model, again consisting of the lithosphere, the upper mantle, and the upper and lower parts of the lower mantle, with the interfaces at depths of 220, 670, 1500 and 2890 km. Now the conductivities of the UM and ULM are considered to be ϑ-dependent, such that the cut-off degree J in the conductivity parameterization (54) is equal to J = 1. The conductivity of the lithosphere is again fixed to 0.001 S/m. Because of the rather poor resolution of the LLM conductivity, as demonstrated in the previous section, this conductivity is chosen to be equal to the optimal value obtained by the CG minimization, that is 97.8 S/m, and is kept fixed throughout the sensitivity tests and subsequent inversion. Complementary to the sensitivity tests for the zonal coefficients σ 0 shown in Fig. 7, we now carry out the sensitivity tests for non-zonal coefficients σ 1 of the ULM ( = 1) and UM ( = 2) in a way similar to that applied in Section 6.1.1 with the same nominal values for the zonal coefficients σ 0 and σ 0 11 = σ 0 21 = 0. The forward and adjoint solutions are again computed for the 2001-CHAMP data (see Section 4) with spherical-harmonic cut-off degree j max = 4 and time step t = 1 h. The earth model is again divided into 40 finite-element layers with layer thickness increasing with depth.
In Fig. 9, we summarize the results of sensitivity tests. The top panels show the misfit χ 2 as a function of the parameters σ 11 and σ 21 , where only one conductivity parameter is varied and the other to zero. The bottom panels compare the derivatives of the misfit obtained by the ASM and the BFM. We can see that the adjoint sensitivities show very good agreement with the brute-force results, with differences not exceeding 0.01 per cent of the magnitude of the sensitivities themselves. Moreover, the sensitivities to latitudinal dependency of conductivity are significant, again more pronounced in the upper mantle than in the lower mantle. This tells us that the CHAMP data are capable of revealing lateral variations of conductivity in the upper and lower mantle.

Conjugate gradient inversion
The sensitivity results in Fig. 9 encourage us in attempting to solve the inverse problem for lateral variations of conductivity in the mantle. For this purpose, we again employ the CG minimization with derivatives obtained by the ASM. The inverse problem is solved for four parameters, σ 0 and σ 1 , = 1, 2. The starting values of σ 0 are the nominal values of the three-layer, 1-D conductivity model (see Section 5.3), while the values of σ 1 are put equal to zero at the start of minimization.
The results of the inversion are summarized in Fig. 10, where the left and centre panels show the conductivity structure in the ULM and UM, while the right panel shows the misfit χ 2 as a function of CG iterations. The blue lines show the starting model of minimization, the dotted lines the model of minimization after the first iteration and the red lines the final model of minimization after 8 iterations. These models are compared with the optimal three-layer, 1-D conductivity model (black lines) found in Section 6.1.2 Again, as indicated by the sensitivity tests, the minimization, at the first stage, adjusts the conductivity in the upper mantle, to which the misfit χ 2 is the most sensitive, and then varies the ULM conductivity, to which the misfit is less sensitive. The optimal values of the conductivity parameters after eight iterations are (σ 10 , σ 11 , σ 20 , σ 21 ) = (0.192, −0.008, −0.476, 0.106). We conclude that the mantle conductivity variations in the latitudinal direction reach about 20 per cent of the mean value in the upper mantle and about 4 per cent in the upper part of the lower mantle. Comparing the optimal values of the zonal coefficients σ 10 and σ 20 with those found in Section 6.1.2 for a 1-D conductivity model, we also conclude 1388 Z. Martinec and J. Velímský  that the averaged optimal 2-D conductivity structure closely approaches the optimal 1-D structure. This is also indicated in Fig. 10, where the final 2-D conductivity profile (red lines) intersects the optimal 1-D conductivity profile (black lines) at the magnetic equator.

Low-pass spatial and frequency filtering of CHAMP data
So far, the sensitivity analysis has been performed for the time series of CHAMP spherical harmonic coefficients X (obs) j (t) and Z (obs) j (t) up to spherical degree j max = 4. However, the traditional approach of inferring the mantle conductivity structure from magnetic time series considers only the first-degree spherical harmonics, X (obs) 1 (t) and Z (obs) 1 (t), or their linear combinations (e.g. Kuvshinov & Olsen 2006;Velímský et al. 2006). We are now interested in the question of how the cut-off degree of the spherical harmonics influences the sensitivity of the misfit with respect to the mantle conductivity structure.
In order to assess this effect, we reduce the cut-off degree of the spherical harmonic expansion of the CHAMP time series to j max = 1. The resulting adjoint sensitivities of misfit χ 2 for the three-layer, 1-D mantle conductivity model employed in Section 5.3, are depicted in Fig. 11. Comparing it with Fig. 7 for j max = 4, we can see that decreasing the cut-off degree of the CHAMP data, which means applying a low-pass spatial filter, reduces the sensitivity of the misfit with respect to the upper mantle conductivity, whereas the sensitivities with respect to lower-mantle conductivities are changed only slightly. This is explained by the fact that the spherical harmonics of the Green's function of EM induction for an external excitation become localized closer to the Earth's surface when increasing their spherical degree. Fig. 11 also shows that the misfit for j max = 1 is mostly sensitive to conductivity variations in the upper part of the lower mantle, which was also demonstrated by Velímský et al. (2006). Figure 11. As for Fig. 7, but for low-pass spatially filtered CHAMP data with the cut-off degree j max = 1.
In the context of the approximation error of the ASM, a question arises as to how low-pass frequency filtering of CHAMP time series influences the sensitivities of the misfit χ 2 . To assess the effect of low-pass filtering of CHAMP time series on sensitivities, we apply the Nutall four-term cosine low-pass filter (Marple 1987) with a cut-off period of 1.9 days on the CHAMP time series of X (obs) j (t) and Z (obs) j (t) for j = 1, · · ·, 4. Fig. 12 shows the unfiltered and low-pass filtered Fourier amplitude spectrum of X (obs) 1 (t) (red lines) and Z (obs) 1 (t) (black lines), respectively. Similarly to the maximum-entropy power spectrum (Fig. 3), we can see spectral peaks at higher harmonics of the 27-day solar rotation period. The low-pass filtered CHAMP time series of X (obs) j (t) and Z (obs) j (t) are then used to compute the adjoint sensitivities for the 3-layer, 1-D conductivity model, studied in Section 5.3. Comparing the resulting sensitivities depicted in Fig. 13 with those for the CHAMP unfiltered time series shown in Fig. 7 we can see that applying the low-pass filter significantly reduced the approximation error O(∂ 2 A/∂t 2 t) of the adjoint sensitivities with respect to the upper-mantle conductivity, because of the reduction of the magnitude of ∂ 2 A/∂t 2 . There is, however, a price to be paid for reducing this error, namely, the sensitivity of the misfit χ 2 with respect to the upper-mantle conductivity is also reduced. On the other hand, the sensitivities with respect to lower-mantle conductivities are hardly changed. This is explained by the skin effect, which, for a given conductivity structure, is the fast that the depth to which an electromagnetic signal can penetrate decreases as the oscillation frequency of the signal increases. Hence, removing the high-frequency part of the CHAMP time series by low-pass filtering reduces the 'illumination' of the upper mantle by an electromagnetic signal and, thus, reduces the sensitivity of EM induction data to upper-mantle conductivity variations. Velímský et al. (2006) interpreted the CHAMP magnetic data recorded during eleven storm events in the time interval 2001-2003 in terms of Earth's 1-D electrical conductivity structure. They calculated the partial derivatives of the misfit with respect to conductivity parameters by the BFM. Although only eleven shorter time intervals of the CHAMP time series were interpreted, the inverse modelling required long computational time.

C O N C L U S I O N S A N D F U T U R E W O R K
This paper has been motivated by efforts to find an advanced technique for interpreting the time series of CHAMP magnetic data such that the complete time series, not only their parts, can be considered in inverse modelling and still be computationally feasible. It turned out that these criteria are satisfied by a highly efficient method of sensitivity analysis-the adjoint method. We demonstrated this for the year 2001 CHAMP time series with a time step of 1 hr. To apply the adjoint sensitivity analysis to longer time series is straightforward, leading to memory and computational time requirements that are linear with respect to the number of time steps undertaken. The analysis of the complete, more than 8-year long, CHAMP time series is ongoing with the particular objective of detecting the perovskite/postperovskite conductivity contrast in the boundary layer between the lower mantle and outer core. The results of this effort will be the subject of a forthcoming publication. Velímský et al. (2006) showed that the CHAMP data for the period of magnetic storms are most sensitive to conductivity changes in the upper part of the lower mantle. This result was explained by missing higher-frequency signals in the chosen eleven CHAMP records. The sensitivity analysis shown here in Fig. 13 confirms this explanation. Moreover, if a time series longer than the record of a magnetic storm can be taken into account, the spherical harmonic analysis of CHAMP data can be carried out up to higher spherical degrees than degree j = 1 as considered by Velímský et al. (2006), for example, up to degree j = 4 for the 2001-CHAMP data. The sensitivity analysis for the complete 1394 Z. Martinec and J. Velímský In particular, for a toroidal vector A (ϑ), the coefficients A j±1 j (r ) = 0 and eq. (A15) reduces to

A P P E N D I X B : W E A K F O R M U L AT I O N O F T H E A D J O I N T M E T H O D
In this section, we will demonstrate that the adjoint method derived in a strong (differential) form (see Section 3.6) can be reformulated in a weak sense. Such a generalized formulation have all advantages of a weak formulation of any boundary-value problem for the secondorder differential equation (e.g. Křížek & Neittaanmäki 1990), and is particularly suitable for a numerical implementation of finite-element technique which is also applied in Section 5. We further aim at showing that the weak formulation can be set up in a unified manner for both the forward and adjoint methods of EM induction for CHAMP magnetic data. Let us start with the parameterization of toroidal vector potential A in the radial direction. Inside a conducting sphere B, of radius a, the spherical harmonic expansion coefficients A j j (r , t) can be parameterized by P + 1 piecewise-linear finite elements ψ k (r ) on the interval 0 ≤ r ≤ a such that In an insulating atmosphere A, the toroidal vector potential A satisfies the Laplace equation, and thus the spherical harmonic expansion coefficients A j j (r , t) can be parameterized in terms of zonal scalar solid spherical harmonics r j Y j (ϑ) and r − j−1 Y j (ϑ) as follows: where G (e) j (t) and G (i) j (t) are the zonal scalar-magnetic Gauss potential coefficients introduced by eq. (2). Note that the same parameterizations as shown by eqs (B1) and (B2) are taken for the coefficientsÂ j j (r, t). The IBVP (3)-(5) for the satellite magnetic data B t has been formulated in a weak sense by Martinec & McCreadie (2004). We will slightly modify this formulation such that the weak formulation of the forward and adjoint problems will have the same final algebraic forms, and thus implementable in the same numerical way.
The modification concerns the way of how to involve the continuity of toroidal vector potential A on the boundary ∂B between the conducting sphere B and the non-conducting atmosphere A. As seen from parameterizations (B1) and (B2), we apply a different parameterization of A in the sphere B and the spherical layer A. That is why the continuity of A on ∂B must be ensured: where A 0 denotes the potential A in the non-conducting atmosphere A. Martinec & McCreadie (2004) implemented this continuity condition in the construction of a solution space for A in the atmosphere. Here, we apply an alternative approach based on the Lagrange multiplier method.
We introduce the solution spaces V and V 0 for the conducting sphere B and the non-conducting atmosphere A, respectively, as follows: V := { A| A ∈ L 2 (B), curl A ∈ L 2 (B), div A = 0 in B} , where L 2 (B) is the space of square-integrable functions in the domain B and C 2 (A) is the space of functions whose derivatives up to the second order are continuous in the domain A. Note that the space V 0 differs from that used by Martinec & McCreadie (2004), such that the continuity condition (B3) is now not imposed on the functions from V 0 . Instead, we introduce the Lagrange multiplier vector λ and a solution space for it: where L 2 (∂B) is the space of square-integrable functions on the boundary ∂B. Following the considerations of Martinec (1997) and Martinec & McCreadie (2004), the weak formulation of the forward IBVP (3)-(5) consists of finding the potentials A ∈ V and A 0 ∈ V 0 and the Lagrange multipliers λ ∈ V λ such that, at a fixed time t ∈ (0, T ), they satisfy the following variational equality: The bilinear forms a(·, ·), b(·, ·), a 0 (·, ·), c(·, ·) and the functional F(·, ·) are defined as follows: