Applying inversion techniques to derive source currents and geoelectric fields for geomagnetically induced current calculations

This research focuses on the inversion of geomagnetic variation field measurement to obtain source currents in the ionosphere. During a geomagnetic disturbance, the ionospheric currents create magnetic field variations that induce geoelectric fields, which drive geomagnetically induced currents (GIC) in power systems. These GIC may disturb the operation of power systems and cause damage to grounded power transformers. The geoelectric fields at any location of interest can be determined from the source currents in the ionosphere through a solution of the forward problem. Line currents running east–west along given surface position are postulated to exist at a certain height above the Earth’s surface. This physical arrangement results in the fields on the ground having the magnetic north and down components, and the electric east component. Ionospheric currents are modelled by inverting Fourier integrals (over the wavenumber) of elementary geomagnetic fields using the Levenberg– Marquardt technique. The output parameters of the inversion model are the current strength, height and surface position of the ionospheric current system. A ground conductivity structure with five layers from Quebec, Canada, based on the Layered-Earth model is used to obtain the complex skin depth at a given angular frequency. This paper presents preliminary and inversion results based on these structures and simulated geomagnetic fields. The results show some interesting features in the frequency domain. Model parameters obtained through inversion are within 2 % of simulated values. This technique has applications for modelling the currents of electrojets at the equator and auroral regions, as well as currents in the magnetosphere.


Introduction
Solar events, such as coronal mass ejections that become geo-effective, create disturbances within the Earth's magnetosphere giving rise to geomagnetic storms and substorms.During geomagnetic storms, the compression of the magnetosphere by the solar wind and the interaction of the solar wind with the Earth's magnetic field enhance the currents in both the magnetosphere and ionosphere.These currents cause fluctuations in the electric and magnetic fields on the ground.The equatorial electrojet is at 100 km in the ionosphere.The ionosphere extends up to ∼ 1000 km and has current systems that lie at a height of ∼ 100 km.Rapid changes in the geomagnetic field generate geoelectric fields that drive geomagnetically induced currents (GIC) in power lines.The GIC have the potential of causing transformers to fail, with subsequent consequences of a power blackout to the general public, who are increasingly reliant on electrical power for their everyday operations and living (Albertson et al., 1993;Shea and Smart, 1996;Wilkinson, 2007;Withbroe, 2001).
Therefore, it is of interest to power utility operators that a warning system be developed that can predict GIC, which may occur after an eruptive event occurs on the sun.Because of the complexities involved in such a solar-terrestrial interaction and the tremendous challenges facing such a project, we consider as a first step the inversion of the geomagnetic field observations to obtain ionospheric source currents.From these source currents, we estimate the induced geoelectric fields as measured at any location of interest, particularly the electric fields responsible for GIC in power grids on the ground.
The geoelectric fields can be determined from geomagnetic data and the surface impedance (Dearholt and McSpadden, 1973, p. 397), but this calculation is valid only at the location of the known surface impedance.Measured geomagnetic data are readily available for many locations on Earth (Kerridge, 2001;INTERMAGNET, 2014) where there are magnetometers, but conductivity structures for the calculation of surface impedances are scarcer.Therefore, geoelectric data obtained by the direct method are valid locally only (see Fig. 1).
The inversion method allows us to compute geomagnetic fields, via Fourier integrals, at any location of interest.For this to work, a 1-D planar conductivity structure obtained for one location is assumed to be valid everywhere over the relevant region (see Fig. 1 again).From the conductivity structure, a surface impedance or skin depth is computed.
The motivation for using the field inversion method described in this paper is that the geomagnetic measurement for calculating the geoelectric field that drives GIC in power systems is generally not available at the location of interest.Once the parameters of the current system are determined by the inversion from geomagnetic data, one can return to the forward problem Fourier integral and use these parameter values to calculate the geomagnetic field anywhere (the curve in Fig. 1).Inversion provides an alternative way in which to estimate the geomagnetic fields where it is not possible by other means.
Either the electric field can be obtained through a direct multiplication of the spectrum of the magnetic field with the relevant components of the surface impedance at the location of interest, or via the Fourier integral, using a reflection coefficient derived from a model of the ground conductivity.It can be shown that the electric fields are the same when the exact expression for the reflection coefficient is used in the Fourier integrals.But, when the reflection coefficient is replaced with an exponential approximation to facilitate the analytic solution of the integrals, the electric fields obtained are actually different.These differences are addressed in this paper.
A general theoretical framework for computing the fields due to an ionospheric electrojet above a layered Earth was proposed by Hakkinen and Pirjola (1986).To simplify computations, we use the complex image method (CIM), introduced by Wait and Spies (1969) and used by Thomson and Weaver (1975) to replace a conducting planar model for layered Earth by an image current placed at a depth equal to the height of the current system above the Earth plus twice the complex-valued skin depth associated with electromagnetic waves penetrating the Earth.Pirjola and Viljanen (1998) took this approach and applied it to a finite auroral electrojet with field-aligned currents carrying away excess charges to or from the magnetosphere in the Northern Hemisphere.
This paper takes the simulated magnetic data from Boteler et al. (2000) and use inversion techniques in the complex image method to obtain the current strength (and position) of the modelled auroral electrojet.For this purpose we take this ionospheric current to be a line current above the Earth.The planar model in the Cartesian coordinate system with a line current above a layered Earth is a local approximation of a ring current around the Earth.No field-aligned currents are considered here.

Theory
We introduce an inversion approach on simulated magnetic data to obtain ionospheric current system characteristics.The application of the CIM allows one to approximate the reflection coefficient to an exponential.It is dependent on the skin depth, and thus the surface impedance.The surface impedance can be computed for one fixed frequency by using the Quebec conductivity structure in Canada (see Hakkinen and Pirjola, 1986), which is based on the general theory for computing the geomagnetic and geoelectric fields due to an electrojet in the magnetosphere above a layered Earth.Appendix A contains details of the relevant derivations of the theory.
We restate the expressions of Eq. (A9) here: where x is the latitude and ω is the frequency.The parameters are I the ionospheric current strength, h the height, x o the reference latitude.The constant of free-space permeability is µ.
The function p (ω) is the complex skin depth that depends on the conductivity structure.The transpose of a vector is indicated by the symbol T .The Fourier integral expression of ground magnetic field components (Eq. 1) is the model function in the inversion problem that obtains the output parameters (I , h, x o ) of the model for a line current system above the Earth.These parameters, once found, are substituted into the ground geoelectric field expression for estimation: (2)

Forward procedure
A model comprises entities and relations defined by variables and parameters.Entities include the currents, the fields and the Earth.The relations between these entities in the forward problem are Eqs.( 1) and ( 2).The input variables to these equations are the surface distance x and frequency ω.
Output variables are the electric and magnetic fields.The skin depth p is dependent on frequency ω and two sets of parameters: conductivities σ n plus half-space conductivity σ N +1 , and thicknesses h n for level n = 1, . .., N of a layered structure of the Earth.The other parameters are the current strength I and height h.The equations have been altered to include another parameter: the position x o of the current system.Then, x is replaced by x − x o , here.Therefore, the parameter x o shifts the fields along the surface in either direction.
Before any inversion can be performed, reference magnetic data must be obtained against which the inversion can be tested.The test data were obtained by calculating Eq. (1) to replicate the physical set-up of Boteler et al. (2000).In this reference, the physical set-up was a Cauchy distributed current system 100 km above the Earth with a spread of 200 km.Thus, the line current system for this study should be 300 km high to produce the same results.These data have a surface range from −1000 to +1000 km with a grid spacing of 50 km.They can be regarded as a string of magnetometer stations at positions x i along a meridian.Take note that the modelled curves are symmetric around x o = 0 km for B x and E y , and antisymmetric for B z .The current strength was assumed to be 10 3 kA.
For a Quebec conductivity structure and fluctuation period of τ = 5 min, skin depth p(τ ; [h n , σ n ]) and impedance Z(τ ; [h n , σ n ]) estimates were obtained and passed along to the field Eq. ( 1).Here the skin depth and impedance depend on the following variables: τ is the fluctuation period and [h n , σ n ] is the set of conductivities σ n plus half-space conductivity σ N +1 , and thicknesses h n from the conductivity structure.
The results of the forward computations are given in Sect. 4.

Inversion fundamentals
Usually we have a data set d and a model design set m related to each other by an operation F through the relation d = F (m).This defines the forward problem.We only have available the observations d.The process has to be inverted for m = F −1 (d) = G (d), and that requires optimisation techniques and an objective function G(d).This defines the inverse problem.For a more comprehensive description of the theory, consult Chave and Jones (2012).See also Taranatola (2005).
The linear inverse problems take the form d = Fm or m = Gd in which case F and G are matrices and m and d are column vectors.These matrices are constant with respect to m and d.However, the optimisation is in general of a nonlinear nature.Some approaches take advantage of linear methods by considering a linearised form of the inverse problem.This is accomplished by expanding F (m) in a Taylor series around a reference model m * : where J m * is a Jacobian matrix with J m * ij = ∂F i (m) /∂m j m=m * .The • is the norm of a vector and o is the Landau operator from asymptotic theory on the norm of the model difference m − m * .A linearised inverse problem F is formed when only the first two terms in Eq. (3) are retained and the higher-order terms are discarded.F is then an affine transformation: a linear transformation plus a constant.

Errors and standard deviations
Adding an error vector to the data gives d = F (m)+ .This alters the inversion relation to include errors in the model: m = G d .The notion of well-posed problems (forward and inverse) was established by Hadamard (1902).The conditions for well-posed are a solution to a problem must (1) exist, (2) be unique, and (3) be stable.Failure of any one of these conditions results in an ill-posed problem.Thus, we have a forward problem d = F m .The F maps a subset of vectors in model space (the domain of F) to a subset of vectors in the data space (the range of F).The existence of m means that d must be in the range of F. The uniqueness of m follows when F has a one-to-one transformation, mapping different vectors in model space to different vectors in data space.Then the solution is given by m = F −1 d , where F −1 is defined such that its domain is in the range of F. The stability of m pertains to the effect of the error on m.When error-free ( = 0), then m = F −1 (F (m)) = m; that is, the solution is not only unique, it is correct.In general, however, Stability requires that the solution error be bounded when the data error is bounded as well.Thus, m is stable when a positive function ε(µ) exists, such that δm < µ whenever < ε(µ).This is a definition for the continuity transformation: m is stable when F −1 is continuous.
It is possible to calculate variances and standard deviations of the output parameters in the model.The variances and deviations are obtained from a covariance matrix for parameters.We start with the error vector of residuals = [r 1 , r 2 , . .., r N ].The residuals r i are computed by taking the difference between the data and the forward problem function, where N is the number of data.A Jacobian matrix J is formed by partial differentiating of the objective function with respect to each of its parameters for a set of data points x i .Then the deviation of the fitted m from the actual m parameter vector position for the minimum of the objective function is Neglecting higher power terms of , we multiply µ with its transpose: The sum-of-squared residuals (SSR) are obtained from the corresponding multiplication of with its transpose: Here the SSR variation is s SSR .Substituting this into Eq.( 6) and dividing by N, the parameter covariance matrix is formed: The variance for the parameters can then be obtained from the diagonal elements of .The square roots of the diagonal elements are the parameter standard deviations.

An optimisation problem
An optimisation in a simple case is a minimisation or maximisation of a function describing some system characteristic (say a physical property) dependent on m.In an advanced case the objective function f(d, m) might then be subject to equality constraints and/or parameter bounds m L and m H .A general problem description may be stated as follows: This is a minimisation problem.Most optimisation techniques are designed to be minimisation techniques.Maximising the objective function instead requires that function to be negated and the negative function be then minimised again, e.g.f(d, m) = [−f (d, m)].The inequality constraints may be negated as well, that is An efficient and accurate solution to this problem depends not only on the size of the problem in terms of the number of constraints and model design parameters, but also on the characteristics of the objective function and constraints.When both f(d, m) and f η,µ (d, m) are linear functions of the model vector, the problem is known as a Linear Optimisation (LO) problem.Quadratic Optimisation (QO) concerns the optimisation of a quadratic objective function with linear constraints.For both types of problems, reliable solution procedures are readily available, such as the decomposition methods.More difficult to solve are Nonlinear Optimisation (NO) problems, in which f(d, m) and f η,µ (d, m) can be nonlinear functions of the model vector.A solution of the NO problem generally requires an iterative technique to establish a search direction.This is usually achieved by an approximate solution of an LO, a QO or unconstrained sub-problem.

Least-squares problems
When the optimisation problem is a least-squares problem, the objective function f(d, m) assumes the form of a sum-of-squares function of residuals.That is for data d The same constraints apply and f η,µ (m) are arbitrary functions.Linear least squares will not be used in this study, so decomposition techniques does not apply.One has to rely on iterative optimisation algorithms.Many iterative techniques can be applied on nonlinear least-squares inversion problems.These require much computational work, representing the different methods in which a nonlinear model starts at an initial guess position m s , and is brought closer to the position m m of a minimum of the objective function by an appropriately determined search vector s at each iteration.The m s can be arbitrarily chosen by the user, but it should be in the neighbourhood of a local m o to ensure convergence of that minimum.Otherwise, the technique converges to a wrong minimum or does not converge at all.When convergence is too slow, the technique stops after a maximum number of iterations has been reached and then outputs a warning.

Levenberg-Marquardt algorithm
The Levenberg-Marquardt (LM) algorithm is a restrictedstep method, in only the L 2 norm for least-square nonlinear problems that locates a minimum of a function expressed as the sum of squares of nonlinear functions.According to the abstract of Lourakis ( 2005) "[It] can be thought of as a combination of [the] steepest-descent and the Gauss-Newton method[s]".When the current solution is far from the correct one, the steepest-descent behaviour dominates: slow but guaranteed to converge.When close to the correct solution, the Gauss-Newton behaviour takes over.
We map an output parameter vector m ∈ R m to a measurement vector d ∈ R n with an assumed function d = F(m).An initial parameter estimate m 0 and corresponding measurement x is provided.It is desired to find a vector m min that best satisfies the functional relation F, i.e. that minimises the squared distance T , with = d − d = δd.The basis of LM is a linear approximation to F in the neighbourhood of m.The symbol || • || is a 2-norm.For a small ||δm|| a Taylor series expansion leads to the approximation where J is the Jacobian matrix of F(m).Like all nonlinear optimisation methods, LM is iterative -starting from p 0 it produces a series of vectors m 1 , m 2 , . . .that converge to the local minimiser m min for F. At each step, it is required to find the model change δm that minimises The desired δm is therefore a solution to a linear least-square problem: the minimum is obtained when Jδm − is orthogonal to the column space of J. Thus, J T (Jδm − ) = 0.This yields δm as a solution to the so-called normal equation: The matrix J T J on the left-hand side of Eq. ( 12) is an approximate Hessian.The LM actually solves a variation of Eq. ( 12), known as the augmented normal equations: where N = J T J + αdiag(J T J).The strategy of adjusting diagonal elements of N is called damping and the α is referred to as the trust-region damping term.
If the updated parameter term leads to a reduction in the error , the update is accepted and the process repeats with a decreased value of α.Otherwise α is increased, the augmented normal equation is solved again, and the process iterates until a value of δm is found that decreases error.The process of repeatedly solving Eq. ( 13) for different values of the damping term until an acceptable parameter vector update is found corresponds to an iteration of the LM algorithm.
If the damping term is set to a larger value, the matrix N is nearly diagonal and the LM update step δm is near the steepest-descent direction.The magnitude of δm is reduced contributing to its slowness in this behaviour.Damping also handles situations where the Jacobian is rank deficient and J T J is singular.The LM then defensively navigates a region of the parameter space where the model is nonlinear.If damping is small, the LM step approximates the exact quadratic step appropriate for a linear problem in a Gauss-Newton way.LM is adaptive because it controls its own damping.It raises damping if a step fails to reduce error.Otherwise, damping is reduced.In this way, the LM is capable of alternating between a slow descent approach when far from the minimum and a fast convergence when in the neighbourhood of the minimum.

The computer software
The inversion set-up used in this study is an optimisation curve-fitting tool in the Matlab programming language (MathWorks Inc., 2012) and the inversion problem is summarised in Table 1 and Fig. 2.
The Levenberg-Marquardt algorithm in Matlab terminates when at least one of the following conditions (exitflag value) is met: 4 Magnitude of search direction is smaller than the specified tolerance.
3 Change in the residual was less than the specified tolerance.
2 Change in m was less than the specified tolerance.
1 Function converged to a solution m.
0 Number of iterations exceeded option "MaxIter", or number of function evaluations exceeded option "Max-FunEvals".
−1 Output function terminated the algorithm.
−2 Problem is infeasible: the bounds are inconsistent.
−4 Optimisation could not make further progress.

Inversion procedure
One needs a space domain over which the inversion must run.That is provided by the surface position x, with a grid of data points d i = d B x , d B z T (x i ), 50 km apart from −1000 to +1000 km, where x i is the ith position of the datum d i along the meridian.The frequency ω can also be a domain over which a different inversion problem could run.However, for purposes of this paper the period of the frequency is fixed to τ = 2π/ω = 5 min.An objective function is the sum-of-squared residuals: The model function is the ground geomagnetic field component expressions of Eq. ( 1).The components are complex valued and the code of the curve-fitting toolkit cannot operate on complex values.Therefore, the model function should be expanded from the usual vector form B x , B z to a matrix ReB x , ReB z ; ImB x , ImB z instead.Here commas separate columns and semi-colons separate rows.The geoelectric field on the ground, Eq. ( 2), remains in the forward problem.Once the set of output parameters are found by inversion, they can be substituted into Eq.( 2) to estimate the geoelectric field.
Using the LM technique the objective function (sum-ofsquared residuals) is minimised to determine the output parameters of the model.In the LM, there are no equality or inequality restrictions, but bounds can be set for the parameters.
The aim of the inversion is to optimise f (d, m) to an input data set d of magnetic values reproduced here in the forward problem.Outputs are the parameter set m from any elements in the current system set [I, h, x 0 ] and the layered Earth set [h n , σ n ] plus σ N+1 .In a full inversion, all the parameters are adjusted simultaneously; otherwise, the inversion is partial with one or more parameters fixed and at least one parameter adjusted.For instance, adjusting only the current, when the other parameters do not take part, is a partial inversion.On the other hand, a full inversion adjusts all the parameters of both sets combined.Depending on the aims and scope of any geophysical research project that involves inversion theory, any combination of any number of parameters from any set can be used in the optimisation (such as m = [I, h n , σ m ]; mn = 1, 2, . ..N).In this study, however, we will concern ourselves only with adjusting the current system set of parameters and fix the layered Earth set to the values of Quebec's structure.Thus, From the placement of the current strength I in the magnetic field equations (Eq.1), it is clear that the current strength is a linear model parameter, leading to a linear least-square inversion problem when only this parameter is adjusted.The current is unbounded and can even be zero or negative.Placement of the height h and surface position x 0 in those field expressions turns the inversion problem into a nonlinear least-square fit.The surface position is also unbounded in both the negative and positive directions.
The height can have no negative values however, hence the lower bound of 0 km (i.e. the surface).The skin depth and impedance are not output parameters to the inversion, as they are dependent on output parameters from the structure set.
Since the structure set is fixed, these two surface quantities will be fixed when the period of τ = 5 min is fixed.Optimisation results are shown in Sect. 5.

Preliminary work
The 1-D approximation of the ground conductivity structure of Quebec, Canada, based on magneto-telluric measurements, is summarised in Table 2. Quebec appears to have a resistive structure.This can be used to calculate the material properties (skin depth and impedance) at the surface.This structure determines how the magnetic and electric fields behave.
The skin depth value at period τ = 5 min is 135.122-80.950ikm (or in terms of complex amplitude and phase, 157.52 km at −30.93 • ) and the impedance value is 2.131 + 3.556i m (or 4.15 m at 59.07 • ).Here the impedance is 90 • ahead of the skin depth.

Geomagnetic and geoelectric fields
Once the impedance and skin depth were evaluated at the given period, one works out the respective electric and magnetic fields (still in the forward problem, and shown in Fig. 3) of a line current with strength 1000 kA, positioned at x o = 0 km and a height of 300 km above the surface of the Earth.The extreme values obtained by reading off from the plots of Fig. 3 are listed in Table 3.
Thus, magnetic component B x oscillates almost in phase with fluctuations in the current, while component B z is almost out of phase with the current (between 29.74 and 8.75 • short of 180 • ).The electric component E y is more than 90 • behind the current (between 51.84 and 78.69 • ahead of 180 • ).
Figure 3 gives a general idea of how the fields behave in the surface position space.These can be used in an inverse problem, for example to narrow down the region of interest and provide reasonable starting points for the search of the optimal point in parameter space.

Inversion results
Using the data reproduced in Fig. 3, an inversion was performed as a test to determine the parameters for the current system.This was done to make sure the inversion works properly and to check that the output parameters settle close to the expected values.The inversion worked no matter how far the parameters were initialised from their expected values (as given in the caption of Fig. 3).The results of a full inversion are given in Fig. 4.
When all the parameters of a model are estimated in the inversion, that inversion is called a full inversion.When some parameters are fixed, that inversion is called a partial inversion.For partial inversions with either or both distance parameters fixed, the distance parameters have the constant values given in the caption of Fig. 3.The parameter of the current was never fixed in all inversion cases.The fitted parameters were initialised to the values given in the caption of Fig. 4.
Table 4 shows the final parameter values after the inversion in the three cases where the current and one or both distance parameters were varied.The full inversion is "Case 1", while the partial inversions are "Case 2" and "Case 3" respectively (with only one fixed parameter).All parameters are within 2 % below their values given in the caption of Fig. 3. Rerunning the inversion with both distance parameters fixed, thus varying only the current, produces the current strength at I = 980±2.405kA (or I /I = 0.245 %).This is not shown in Table 4; but it may be labelled as "Case 4".The residuals are not randomly distributed, as could be expected from a Gaussian distribution of errors.The inversion is nevertheless a close to optimal fit of the model to the data.Inversion output parameter standard deviations, denoted by m in Table 4, are derived in Sect.3.For further information we refer to Chave and Jones (2012).The standard deviations are increasing when inverting Case 1 through Case 4 in that order.Table 4 shows how they increase.A higher deviation means the parameter is more unstable.The fewer the number of current system parameters involved in the inversion, it seems the more unstable they become.The legends are defined as follows: "Data" is simulated magnetic field observations, "Init" is the initial estimate of the magnetic field before the inversion, "Final" is the final estimate of the magnetic field obtained through the inversion.The legends are defined as follows: "Data" is electric field simulations from simulated magnetic field observations via E y (x i , ω) = −Z(ω)H x (x i , ω) (see Fig. 1), "Init" is the initial estimate of the electric field before the inversion (using parameter initial values), "Final" is the final estimate of the electric field obtained after the inversion (using parameter final values obtained by the inversion).

Conclusions
This paper demonstrates the use of inversion techniques, using the complex image method to determine the parameters (current strength and/or two positions) of the line current by fitting the ionospheric currents to magnetic data calculated from Eq. (1).For this purpose the ionospheric current was taken as a line current above the Earth.No field-aligned currents were considered.
The current is the most important of the three parameters involved and always varies in each inversion case.Case 3 shows the best current value (0.2 % error) against a target of 1000 kA.Cases 1 and 2 are intermediate; the current is 1 % below target.Case 4 is the worst case: the current is 2 % below target.
The value of the estimated current height decreases by 35 km when the surface position is fixed (in both Cases 1 and 2) and gets closer to a target value of 300 km.The value of the estimated surface position moves south to within 2 km of the target at the equator, when the height is fixed, but not crossing the equator (in Cases 1 and 3 respectively).
After the inversion from magnetic data, the electric field were inferred as shown in Fig. 5.The estimated electric field can then be used to determine the GIC in power networks.
Variations in the horizontal magnetic field components B x and B y induce a geoelectric field which then drives an electric current in the Earth according to Ohm's law J = σ E. The geoelectric field at the Earth's surface can be modelled using the plane wave model (Viljanen and Pirjola, 1989;Pirjola, 2002).

A2 Homogeneous Earth model
As a first approximation, we assume the Earth is a uniform half-space of homogeneous conductivity and assume that there is a plane wave field that propagates vertically downwards.Using a Cartesian coordinate system where the xy plane corresponds to the Earth's surface, then at a single frequency ω the fields of a plane wave can be expressed as where For the given frequency ω, the propagation constant k is given by k = −iκ, where κ = ω √ ε o µ o .For a lossy medium the skin depth is complex.For a good conductor and low frequencies the quasi-static approximation (ωε o /σ 1) can be applied: where µ o and ε o are the permeability and permittivity constants of free space and σ the uniform conductivity.Here use was made of a complex identity In a homogeneous conducting medium with uniform conductivity, the plane wave amplitude decays with depth into the medium.The depth at which the amplitude has decayed to e −1 times the amplitude at the surface is the skin depth δ = √ 2/ωµ o σ .And the complex skin depth is p.It can be shown that κδ = √ 2ωε o /σ .The approximation of the reflection coefficient by an exponential function is based on the assumption κδ 1.This assumption is justified in the context of GIC modelling since the spectrum of the geomagnetic field is typically in the range 1 to 10 mHz.For a homogeneous ground conductivity of σ = 1 mS m −1 , which is typical for the locations of interest, the values of δκ are in the range 10 −6 to 10 −5 .
With Eq. (A2) substituted into Eq.(A1) it follows that the ratios between the orthogonal electric and magnetic field components define the surface impedance Z(ω).The ratio between the electric field and the spectral component at ω of the time derivative of the magnetic field, ∂B ∂t , which will be denoted by Ḃ (referred to as B-dot) defines the complex skin depth p (Deri et al., 1981).Thus,

A3 Elementary fields
From Maxwell's equations in the plane-Earth model (Fig. A1) and the quasi-static approximation, a diffusion equation is derived and an electric field elementary solution is found (Hermance and Peltier, 1970) There are however both incident e −γ z and reflected e +γ z waves; the solution is symmetrical in x around x = 0, and above the Earth's surface (z > 0)γ = ν because σ = 0 there.Thus, the electric field is given by where R(ω) is the reflection coefficient and C is an arbitrary constant.
The magnetic field components are computed by taking the curl of the diffusion equation and using Maxwell's equations.It then follows that the only nonzero components are

A4 Layered-Earth model
We next consider a multi-layered model of the Earth.The appendix of Wait (1980) described a general approach to determine Z(ω) from the 1-D multilayer ground conductivity structure of a given location.We assume N planar layers in the ground below the Earth's surface.Each layer (n = N, . .., 1) has a finite thickness h n and a uniform conductivity σ n .Correspondingly, uniform elementary impedances can be obtained from the conductivity for each layer.We define a modified wave number for each layer κ n .Equation (A3) still applies, but the µ o , ε o and σ are replaced by µ n , ε n and σ n respectively for each layer.The intrinsic layer impedance is defined by K n = iωµ /κ n and related to the layer reflection by R n = K n −Z n K n +Z n .For a good conductor in quasi-static approximation (ωε n /σ n 1), we have κ n = √ iωµ n σ n and K n = √ iωµ n /σ n .The (N + 1)th layer is called the remaining half-space in plane-Earth geometry and is assumed to have infinite thickness, uniform conductivity σ N +1 and layer impedance K N +1 .These layer impedances are independent of each other.To relate them, a second set of impedances Z n at the boundaries are defined which are dependent on the layer thicknesses and layer impedances of the layers below and up to that boundary.For the lowest boundary, separating the half-space from the next layer, the boundary impedance is Z N+1 = K N +1 .Thus, a recursion relation is set up, starting at the bottom and working all the way up to the top (that is for n = N, . .., 1): Then the surface impedance is the boundary impedance of the Earth's surface: Z(ω) = Z 1 .In general, the constants of permittivity ε n and permeability µ n are all different for each layer.In the present study the layer permittivities are all set to ε n = ε 0 and the layer permeabilities to µ n = µ 0 for all n.

A5 Complex image method
Next, we consider the complex image method (CIM) and an approximation to the reflection coefficient to accommodate this method.For convenience, this also introduces an equally important material property called the skin depth for multi-layered Earth.The other important material property is the surface impedance.The surface skin depth is computed from the surface impedance as p(ω) = Z(ω)/iωµ 0 .However, the surface is a boundary of the layered-Earth model and, as with boundary impedances Z n , one can also form a set of boundary skin depths p n similarly computed from these impedances (with µ n replacing µ 0 ).
The reflection coefficient (Boteler and Pirjola, 1998) can be expressed as Note that R depends not only on angular frequency ω and wave number ν, but also on the complex surface impedance Z(ω) or skin depth p(ω).
Under the condition that pν 1, it can be shown that the reflection coefficient can be written in exponential form which facilitates analytic solution of the inversion integrals.Replace the Taylor expansion of Eq. (A7) with the Taylor expansion of an exponential function; then, R ≈ e −2pν .This can then be inserted into the Fourier integral expressions for the magnetic and electric fields, which then makes it possible to derive their solutions analytically.
Here the image current is employed to represent the reflected part of the electromagnetic field off the Earth's surface (or equivalently a layered conductive Earth).An image line current is assumed to be flowing in the opposite direction to the external line current at a depth z = h + 2p.

A6 Geomagnetic and geoelectric fields in a plane-Earth model
To relate the elementary fields (in Sect.A3) to that of a line current, one must take Fourier integrals of the components over propagation space ν.This forms the total fields over surface distance and frequency space at the Earth's surface (z = 0).Adapted from Boteler et al. (2000), the geoelectric and geomagnetic field components are then  Boteler et al. (2000) discuss distributions of currents of one type and points out a field equivalence.For a current system defined by a Cauchy distribution, characterised by a spread parameter a, we have j (x) = I π a a 2 +x 2 and the distribution of currents in propagation space is J (ν) = I e −|ν|a .When J (ν) is replaced in the integrals, the extra exponential factor produced is absorbed into the exponential factor of the height: e −νa e −νh = e −ν(h+a) .Thus, a Cauchy spread parameter a is then added to the height h.An equivalence of representation has resulted: the fields produced by using a Cauchy distributed current placed at a height z = −h, would be equivalent to those created by a line current system placed at a height of z = − (h + a).The new height can be denoted as z = −h .Cauchy distributed current systems are represented by line currents further from the z = 0 interface (the Earth's surface), as determined by the Cauchy parameter.Therefore, one can disregard the need for such distributions, and consider only line current systems.A line current would only have j (x) = I δ (x), leading to J (ν) = I ∞ −∞ δ (x) e −iνx dx = I .The integrals will need to be solved numerically if the exact expression for R in Eq. ( A7) is substituted for the reflection coefficient.No closed analytic solutions exist for the combination of elementary functions present in the resulting integrands.However, replacing the reflection coefficient in the integrals by its approximation to Eq. (A7) means there will only be two elementary functions in the integrands: the trigonometric and exponential functions.Exact solutions for these types of integrals have been derived, and that serves as a motivation for using the image current method in simplifying the derivation and evaluation of these integrals.Making all the substitutions to Eq. (A8), the final form solutions can be obtained from any standard integral table, and is given as

Figure 1 .
Figure1.A conceptual diagram to illustrate the following: (1) the surface impedance for a layered medium independent of position x, (2) the electric field E y can be determined from the direct relation E y (x i , ω) = −Z(ω)H x (x i , ω) (shown as stars) based on the measured magnetic field (shown as dots), (3) the electric field and the magnetic field (shown as solid lines) can be determined from the Fourier integrals of the currents obtained from the inversion of the measured magnetic field, (4) the direct and indirect methods give approximately the same values for the electric field.
Figure 2. A sketch of the forward (left) and inverse (right) problem.The residual is the difference between modelled and measured magnetic field: r i = d i − G(x i ; m) for the ith data point d i at surface position x i .Sum of squared residuals over surface position is S = i r 2 i for a least-squares inversion.G is the objective function Eq. (1).

Figure 3 .
Figure 3. Simulated magnetic (B x , B z ) and electric (E y ) field component plots against surface position x in the forward problem (parameters: h = 300 km, x o = 0 km, I = 10 3 kA) for the Quebec structure and period τ = 5 min.Each complex part is plotted separately.All plots are symmetric (B x , E y ) or anti-symmetric (B z ) around x = 0 km.

Figure 4 .
Figure 4. Fitted magnetic fields B x , B z from inversion for the Quebec ground conductivity structure at a period of τ = 5 min.Initial values of the parameters were as follows: current strength I = 5 × 10 4 kA, height h = 6000 km, surface position x o = 1000 km.Target parameter values were as follows: I = 10 3 kA, h = 300 km and x o = 0 km.A positive value of x o means a position north of the equator.All fitted parameters came close to their targeted values.The legends are defined as follows: "Data" is simulated magnetic field observations, "Init" is the initial estimate of the magnetic field before the inversion, "Final" is the final estimate of the magnetic field obtained through the inversion.

Figure 5 .
Figure 5.Estimated electric field E y from inversion of the simulated magnetic field data shown in Fig. 4 using the Quebec ground conductivity structure at a period of τ = 5 min.Initial values of the parameters were as follows: current strength I = 5×10 4 kA, height h = 6000 km, surface position x o = 1000 km.Target parameter values were as follows: I = 10 3 kA, h = 300 km and x o = 0 km.A positive value of x o means a position north of the equator.All fitted parameters came close to their target values.The legends are defined as follows: "Data" is electric field simulations from simulated magnetic field observations via E y (x i , ω) = −Z(ω)H x (x i , ω) (see Fig.1), "Init" is the initial estimate of the electric field before the inversion (using parameter initial values), "Final" is the final estimate of the electric field obtained after the inversion (using parameter final values obtained by the inversion).

Figure A1 .
Figure A1.Plane-Earth model of the current image method for a 1-D representation of a conductivity structure.Labels define the concepts of the method.

Table 1 .
Summary of the inversion set-up.

Table 3 .
Extreme values of the magnetic and electric field components at period τ = 5 min obtained by using the ground conductivity structure of Quebec.Amplitude and phase ∼ (1004.04 at 5.14 • ) ∼ (228.04 at 164.74 • ) ∼ (3.23 at 248.20 • )

Table 4 .
Deviations m from the nominal values m of the inverted parameters in the form m/m, from the inputs of the published data.The notation of m ∈ [I, h, x o ] in each case is the parameter of interest.