Time-domain modelling of 3-D Earth's and planetary electromagnetic induction effect in ground and satellite observations
Alexander V. Grayver, Alexey V. Kuvshinov, Dieter Werthmüller
mmanuscript submitted to
JGR: Space Physics
Time-domain modeling of 3-D Earth’s and planetaryelectromagnetic induction effect in ground and satelliteobservations
Alexander V. Grayver , Alexey Kuvshinov , Dieter Werthm¨uller Institute of Geophysics, ETH Zurich, Sonneggstrasse 5, CH-8092 Zurich Faculty of Civil Engineering and Geosciences, TU Delft, Netherlands
Key Points: • Accurate modeling of EM induction effects in ground and satellite observationsvia local and global impulse responses; • Including 3-D EM induction effects improves description of the observed magneticvariations during both quiet and disturbed conditions; • We provide a data constrained model of external currents for the
Swarm era;
Corresponding author: Alexander V. Grayver, [email protected] –1– a r X i v : . [ phy s i c s . g e o - ph ] N ov anuscript submitted to JGR: Space Physics
Abstract
Electric currents induced in conductive planetary interiors by time-varying magnetosphericand ionospheric current systems have a significant effect on electromagnetic (EM) fieldobservations. Complete characterization of EM induction effects is difficult owing to non-linear interactions between the three-dimensional (3-D) electrical structure of a planetand spatial complexity of inducing current systems. We present a general framework fortime-domain modeling of 3-D EM induction effects in heterogeneous conducting plan-ets. Our approach does not assume that the magnetic field is potential, allows for an ar-bitrary distribution of electrical conductivity within a planet, and can deal with spatiallycomplex time-varying current systems. The method is applicable to both data measuredat stationary observation sites and satellite platforms, and enables the calculation of 3-D EM induction effects in near real-time settings.
The effect of electric (telluric) currents induced in subsurface was observed in timeseries of geomagnetic field variations as early as in Schuster (1889), where it was also pro-posed that this effect depends on the electrical conductivity at depth. Subsequent stud-ies have led to the establishment of an entirely new research field that exploits the elec-tromagnetic (EM) induction phenomenon to sound planetary interiors (Price, 1967). Nev-ertheless, present studies focusing on natural current systems, such as in magnetosphereor ionosphere, often neglect the effect of currents induced in the subsurface or treat itby using a variety of simplistic assumptions. However, as model parameterizations havebecome more realistic and accuracy of the geomagnetic measurements has improved, theeffect of induction may no longer be neglected or substantially simplified, creating a needfor efficient methods which can accurately account for it.There exist two principal approaches to account for the induction effects in geo-magnetic data. First, one can separate an observed vector magnetic field into inducing(external) and induced (internal) parts by using the classic Gauss’ method (Gauss, 1877).However, limitations imposed by this method, namely that the magnetic field must bepotential and measured in a region between the inducing and induced currents, eitherrestrict or invalidate its applicability. Besides, including more unknowns in statisticalmodels to constrain the induced part may quickly degrade the quality of the models givennoisy data with limited coverage. An alternative to that is to model the EM inductioneffect due to extraneous currents by invoking the governing Maxwell’s equations.The latter approach has several advantages. Unlike Gauss method, the modelingapproach is applicable regardless of the position of measurements relative to the induc-ing and induced current regions, and remains valid in regions where the field is not po-tential. Additionally, this has a positive effect for the conditioning of statistical modelssince extra unknowns used to describe the induced part can be eliminated.In practice, the complication behind modeling EM induction in geomagnetic ob-servations is twofold. First, one needs to assume a subsurface conductivity model. A num-ber of regional and global conductivity models exist. This study will not focus on howthese models are constructed and whether they represent the subsurface accurately, eventhough inaccurate conductivity models may bias results. It is expected that our knowl-edge about the electrical structure of the subsurface will continuously improve, allow-ing for the construction of more accurate models at different scales (Kelbert, 2020). Sec-ond, even if a distribution of the subsurface conductivity was known, modeling inducedresponse of a 3-D heterogeneous planet remains a computationally demanding problem.Our goal is to develop an efficient time-domain method for calculating the EM in-duction effect of a planet with an arbitrary 3-D conductivity distribution that is suit-able for both ground and satellite observations. –2–anuscript submitted to
JGR: Space Physics
One way to calculate a planet’s EM induction effect is through frequency domain(FD) transfer functions, which describe a planet’s response due to ”elementary” extra-neous currents. Modeling 3-D EM induction effects with transfer functions was previ-ously applied to analyse daily magnetic field variations (Yamazaki & Maute, 2017) inground (Kuvshinov et al., 1999; Koch & Kuvshinov, 2013; Guzavina et al., 2019) andsatellite measurements (Sabaka et al., 2004, 2015; Chulliat et al., 2016; Sabaka et al., 2018).Additionally, it was applied in the analysis of aperiodic geomagnetic variations in groundobservations (Olsen & Kuvshinov, 2004; P¨uthe et al., 2014; Sun et al., 2015; Honkonenet al., 2018; Munch et al., 2020). These studies Fourier transformed data and appliedtransfer functions in frequency domain, followed by inverse Fourier transform in orderto obtain results in time domain. Therefore, all aforementioned studies effectively workedin frequency domain.However, the FD approach based on transfer functions has limitations in many prac-tical scenarios. Among them are applications involving nearly real-time predictions ofinduction effects with constantly augmented time series, such as space weather hazardassessment, or estimation of steering errors in geomagnetic navigation while drilling. Thelimitations of the FD approach are also apparent when working with data from constantlymoving satellites due to spatio-temporal aliasing. To overcome these restrictions, trans-fer functions can be converted into impulse responses and applied to the data directlyin time domain. This approach was adopted by Maus and Weidelt (2004); Olsen et al.(2005); Thomson and Lesur (2007) for modeling EM induction effects in satellite data.However, these works only considered the induction effect due to an external source de-scribed by a single (first zonal) spherical harmonic function and, moreover, assuming a1-D subsurface conductivity distribution. The extension of this concept to general set-tings and presentation of all methodological details constitute the main contribution ofthis study.Here, we calculate time-domain impulse responses of a medium by converting trans-fer functions pre-calculated in frequency domain. We achieve high computational effi-ciency by applying optimal digital linear filters (DLF) (Ghosh, 1970, 1971a) with the laggedconvolution method (Anderson, 1975), which require only a small set of (computation-ally expensive) frequency domain solutions. For this purpose, we design new DLFs us-ing the methodology presented in Werthm¨uller et al. (2019). Alternatively, evaluationof impulse responses of a 3-D medium can be done by means of dedicated time-domaininduction solvers (Vel´ımsk`y et al., 2003; Vel´ımsk`y & Martinec, 2005).The methods developed here are applied to describe induction effect due to iono-spheric and magnetospheric currents in ground and satellite geomagnetic observations.However, the formalism is amenable to observations made around other planets, whereconventional methods may be too restrictive (e.g. Olsen et al., 2010).
Electromagnetic field variations are governed by Maxwell’s equations. In frequencydomain, these equations read 1 µ ∇ × (cid:126)B = σ (cid:126)E + (cid:126)j ext , (1) ∇ × (cid:126)E = − i ω (cid:126)B, (2)where µ is the magnetic permeability of free space; ω angular frequency; (cid:126)j ext ( (cid:126)r, ω ) theextraneous (impressed) electric current density; (cid:126)B ( (cid:126)r, ω ; σ ) , (cid:126)E ( (cid:126)r, ω ; σ ) are magnetic andelectric fields, respectively; σ ( (cid:126)r ) spatial distribution of electrical conductivity; vector (cid:126)r =( r, ϑ, ϕ ) describes a position in the spherical coordinate system with r , ϑ and ϕ being –3–anuscript submitted to JGR: Space Physics distance from the planet’s centre, co-latitude, and longitude, respectively. Note that weneglected displacement currents and adopted the following Fourier convention f ( t ) = 12 π ∞ ˆ −∞ ˜ f ( ω ) e i ωt d ω. (3)We assume that the current density, (cid:126)j ext ( (cid:126)r, ω ), can be represented as a linear com-bination of spatial modes (cid:126)j i ( (cid:126)r ), (cid:126)j ext ( (cid:126)r, ω ) = (cid:88) i (cid:126)j i ( (cid:126)r ) c i ( ω ) , (4)where (cid:126)j i ( (cid:126)r ) can, in practice, include electric dipoles, current loops (Sun & Egbert, 2012)or be a continuous function.By virtue of the linearity of Maxwell’s equations with respect to the (cid:126)j ext ( (cid:126)r, ω ) term,we can expand total (i.e. inducing plus induced) EM field as a linear combination of in-dividual fields (cid:126)B i , (cid:126)E i , (cid:126)B ( (cid:126)r, ω ; σ ) = (cid:88) i (cid:126)B i ( (cid:126)r, ω ; σ ) c i ( ω ) , (5) (cid:126)E ( (cid:126)r, ω ; σ ) = (cid:88) i (cid:126)E i ( (cid:126)r, ω ; σ ) c i ( ω ) . (6)The (cid:126)B i ( (cid:126)r, ω ; σ ) and (cid:126)E i ( (cid:126)r, ω ; σ ) fields are solutions of the equations1 µ ∇ × (cid:126)B i = σ (cid:126)E i + (cid:126)j i , (7) ∇ × (cid:126)E i = − i ω (cid:126)B i , (8)and, following definitions in Appendix A, represent EM transfer functions of a medium.Therefore, a transfer function of a planet at a position (cid:126)r depends on the subsur-face conductivity distribution and frequency of excitation as well as on the spatial ge-ometry of the current density expressed through the (cid:126)j i term. We now elaborate on the form of the current density term (cid:126)j ext . In this study, weassume that electric currents flow within an insulated spherical shell above the ground.This allows us to collapse any current density distribution within the shell into a cur-rent sheet characterized by a stream function (cid:126)j ext ( (cid:126)r, ω ) = − δ ( r − b )ˆ e r × ∇ H Ψ( θ, φ, ω ) , (9)where a is planet’s radius, b = a + h , with h being the altitude of the current sheet, ∇ H f = 1 r ∂f∂θ ˆ e θ + 1 r sin θ ∂f∂φ ˆ e φ , (10)and ˆ e r , ˆ e θ and ˆ e φ are the unit vectors of the spherical coordinate system. Consequently,we can expand the stream function as a linear combination of spatial modes and scalarcoefficients, that is Ψ( θ, φ, ω ) = (cid:88) i Ψ i ( θ, φ ) c i ( ω ) . (11)Using eqs. (4) and (11), we can rewrite eq. (9) as (cid:126)j ext ( (cid:126)r, ω ) = − δ ( r − b ) (cid:88) i [ˆ e r × ∇ H Ψ i ( θ, φ )] c i ( ω ) . (12) –4–anuscript submitted to JGR: Space Physics
The choice of spatial functions Ψ i is generally problem dependent. In this study,we will adopt spherical harmonic (SH) representation. Then, for an external source, astream function can be written as (Schmucker, 1985)Ψ e ( (cid:126)r, ω ) = − aµ (cid:88) ( n,m ) ∈M n + 1 n + 1 (cid:18) ba (cid:19) n ˜ ε mn ( ω ) S mn ( θ, φ ) , (13)where S mn ( θ, φ ) = P | m | n (cos θ ) exp (i mφ ) (14)is a spherical harmonic (SH) function of degree n and order m with P | m | n being Schmidtsemi-normalized associated Legendre polynomials and M is a set of SH functions withcorresponding complex-valued SH coefficients ˜ ε mn ( ω ).This allows us to rewrite eq. (12) as (cid:126)j ext ( (cid:126)r, ω ) = (cid:88) ( n,m ) ∈M (cid:126)j mn ( (cid:126)r )˜ ε mn ( ω ) , (15)with (cid:126)j mn ( (cid:126)r ) = δ ( r − b ) µ n + 1 n + 1 (cid:18) ba (cid:19) n − ˆ e r × ∇ ⊥ S mn ( θ, φ ) , (16)where ∇ ⊥ = r ∇ H . Accordingly, following eqs. (5)-(6), total electric and magnetic fieldsat a position (cid:126)r can be expressed as (cid:126)B ( (cid:126)r, ω ; σ ) = (cid:88) ( n,m ) ∈M (cid:126)B mn ( (cid:126)r, ω ; σ )˜ ε mn ( ω ) , (17) (cid:126)E ( (cid:126)r, ω ; σ ) = (cid:88) ( n,m ) ∈M (cid:126)E mn ( (cid:126)r, ω ; σ )˜ ε mn ( ω ) , (18)where (cid:126)B mn , (cid:126)E mn are magnetic and electric field transfer functions due to the current den-sity distribution as given by eq. (16). In what follows, we will work with the magneticfield only, although some applications in the field of space weather modeling may takeadvantage of eq. (18) to work with electric fields.Note that eqs. (13)-(18) are only valid for a source that is external relative to theobserver. The equivalent derivations for internal sources (such as, for example, ionospherein satellite data) can be carried out by taking (Schmucker, 1985)Ψ i ( (cid:126)r, ω ) = aµ (cid:88) ( n,m ) ∈M n + 1 n (cid:16) ab (cid:17) n +1 ˜ ι mn ( ω ) S mn ( θ, φ ) (19)instead of eq. (13). In this section, we present methods to calculate EM signals induced by an electriccurrent of the form (16) and measured on the ground or in space.
For reasons that we discussed in the introduction, it is often more convenient towork with data in time domain. Therefore, total magnetic field at a location (cid:126)r and time t can be best described by eq. (17) after its transformation to time domain. Eq. (17) –5–anuscript submitted to JGR: Space Physics can be written in time domain as a convolution integral (see Appendix A for more de-tails) (cid:126)B ( (cid:126)r, t ; σ ) = (cid:88) ( n,m ) ∈M + ˆ t −∞ (cid:104) (cid:126)B m ( c ) n ( (cid:126)r, τ ; σ ) q mn ( t − τ ) + (cid:126)B m ( s ) n ( (cid:126)r, τ ; σ ) s mn ( t − τ ) (cid:105) d τ, (20)where M + is a set of SH functions with non-negative orders ( m ≥ q, s inducing SHcoefficients; (cid:126)B m ( c ) n and (cid:126)B m ( s ) n are impulse responses of a medium for the q mn and s mn co-efficients, respectively. They can be defined as (cid:126)B m ( c ) n ( (cid:126)r, t ; σ ) = − π ˆ ∞ Im (cid:34) (cid:126)B mn ( (cid:126)r, ω ; σ ) + (cid:126)B − mn ( (cid:126)r, ω ; σ )2 (cid:35) sin ( ωt ) dω (21)and (cid:126)B m ( s ) n ( (cid:126)r, t ; σ ) = 2 π ˆ ∞ Im (cid:34) (cid:126)B mn ( (cid:126)r, ω ; σ ) − (cid:126)B − mn ( (cid:126)r, ω ; σ )2i (cid:35) sin ( ωt ) dω. (22)The integrals in eqs. (21)-(22) are evaluated by using the digital linear filter method asexplained in Appendix B. For satellite measurements, using local impulse responses becomes impractical sinceit requires calculating eqs. (21)-(22) for every satellite location. Therefore, to describeEM induction effects in satellite data, we resort to different transfer functions, namely Q -responses and Q -matrices, which enable factorization of spatial and temporal effects.We note, however, that while transfer functions in eq. (17) are valid everywhere, Q -responsesand Q -matrices are valid only in regions where the magnetic field is potential.Recall that if a magnetic field at a position (cid:126)r and time t is potential, we have (cid:126)B ( (cid:126)r, t ; σ ) = −∇ (cid:2) V e ( (cid:126)r, t ) + V i ( (cid:126)r, t ; σ ) (cid:3) , (23)where inducing and induced parts of the potential are given by V e ( (cid:126)r, t ) = a N (cid:88) n =1 n (cid:88) m =0 [ q mn ( t ) cos( mφ ) + s mn ( t ) sin( mφ )] (cid:16) ra (cid:17) n P mn (cos θ )= Re (cid:40) a N (cid:88) n =1 m (cid:88) m = − n ε mn ( t ) (cid:16) ra (cid:17) n S mn ( θ, φ ) (cid:41) (24)and V i ( (cid:126)r, t ; σ ) = a K (cid:88) k =1 k (cid:88) l =0 (cid:2) g lk ( t ; σ ) cos( lφ ) + h lk ( t ; σ ) sin( lφ ) (cid:3) (cid:16) ar (cid:17) k +1 P lk (cos θ )= Re (cid:40) a K (cid:88) k =1 k (cid:88) l = − k ι lk ( t ; σ ) (cid:16) ar (cid:17) k +1 S lk ( θ, φ ) (cid:41) , (25)where N, K are some constants that truncate series. Note that we stated the magneticfield potential using both real-valued and complex-valued notations with the followingrelation between the coefficients, ε mn = q mn − i s mn , m > q | m | n +i s | m | n , m < q mn , m = 0 . (26) –6–anuscript submitted to JGR: Space Physics
The relation between induced (internal in our case) coefficients g lk , h lk and ι lk is derivedin an identical way.We can now rewrite the induced magnetic field (25) using transfer functions insteadof induced SH coefficients. Before presenting the general case, we first consider a casewhen planet’s conductivity distribution is assumed to be 1-D, i.e., σ ( (cid:126)r ) ≡ σ ( r ). In thiscase, each coefficient ε mn induces one internal coefficient of the same degree and order (e.g.Price, 1967). Inducing and induced coefficients can be related via a scalar transfer func-tion called Q n -response. In frequency domain, this relation reads˜ ι mn ( ω ; σ ) = ˜ Q n ( ω ; σ )˜ ε mn ( ω ) . (27)Note that Q n is independent of order m (Schmucker, 1985).Following derivations in Appendix A, transforming eq. (27) to time domain andseparating spatial sine and cosine terms leads to a pair of convolution integrals g mn ( t ; σ ) = Q n ∗ q mn = ˆ t −∞ Q n ( t − τ ; σ ) q mn ( τ )d τ , (28) h mn ( t ; σ ) = Q n ∗ s mn = ˆ t −∞ Q n ( t − τ ; σ ) s mn ( τ )d τ . (29)Subsequently, substituting eqs. (28)-(29) in eq. (25) yields internal magnetic po-tential V i ( (cid:126)r, t ; σ ) = V i ( c ) ( (cid:126)r, t ; σ ) + V i ( s ) ( (cid:126)r, t ; σ ) (30)with V i ( c ) ( (cid:126)r, t ; σ ) = a (cid:88) ( n,m ) ∈M + [ Q n ∗ q mn ] cos( mφ ) (cid:16) ar (cid:17) n +1 P mn (cos θ ) , (31) V i ( s ) ( (cid:126)r, t ; σ ) = a (cid:88) ( n,m ) ∈M + [ Q n ∗ s mn ] sin( mφ ) (cid:16) ar (cid:17) n +1 P mn (cos θ ) . (32)Note that in a 1-D case, inducing and induced expansions are identical, hence we used n, m in eqs. (27)-(32) for all SH coefficients.For a general 3-D conductivity distribution, σ ( (cid:126)r ) in a planet, each coefficient ε mn induces infinitely many internal coefficients (Olsen, 1999). The relation between induc-ing and induced coefficients is then described by a set of transfer functions called Q -matrix˜ ι lk ( ω ; σ ) = (cid:88) n,m ˜ Q lmkn ( ω ; σ )˜ ε mn ( ω ) . (33)An element of the Q -matrix is given by (P¨uthe & Kuvshinov, 2014)˜ Q lmkn ( ω ; σ ) = 1( k + 1) (cid:107) S lk (cid:107) ‹ S (1) (cid:2) B mn,r ( (cid:126)r a , ω ; σ ) − B m, ext n,r ( (cid:126)r a ) (cid:3) S l ∗ k ( θ, φ ) sin( θ )d θ d φ, (34)where ∗ denotes complex conjugation, (cid:126)r a = ( a, θ, φ ) is the position vector at the sur-face of a planet, and S (1) the surface of a ball with unit radius. The radial magnetic field B mn,r is (numerically) computed for a given 3-D Earth’s model induced by a unit ampli-tude (˜ ε mn = 1) SH current source described by eq. (16), and B m, ext n,r ( (cid:126)r a ) = − nS mn ( θ, φ ) (35)is the inducing (external) part of the radial magnetic field. –7–anuscript submitted to JGR: Space Physics
In this case, the internal magnetic potential becomes V i ( c ) ( (cid:126)r, t ; σ ) = a (cid:88) ( n,m ) ∈M + (cid:88) k,l (cid:104) Q lm,qgkn ∗ q mn + Q lm,sgkn ∗ s mn (cid:105) cos( kφ ) (cid:16) ar (cid:17) k +1 P lk (cos θ ) , (36) V i ( s ) ( (cid:126)r, t ; σ ) = a (cid:88) ( n,m ) ∈M + (cid:88) k,l (cid:104) Q lm,qhkn ∗ q mn + Q lm,shkn ∗ s mn (cid:105) sin( kφ ) (cid:16) ar (cid:17) k +1 P lk (cos θ ) , (37)where (cid:88) k,l = K (cid:88) k =1 k (cid:88) l =0 . (38)After some algebra, impulse responses in eqs. (36)-(37) can be calculated via sinetransform (see eq. A9) of the spectra, which are related to the frequency domain Q -matrix(eq. 33) via equations below (the dependence on ω and σ is omitted).For l > , m >
0: ˜ Q lm,qgkn = ˜ Q lmkn + ˜ Q l − mkn + ˜ Q − lmkn + ˜ Q − l − mkn , (39)˜ Q lm,qhkn = i ˜ Q lmkn + ˜ Q l − mkn − ˜ Q − lmkn − ˜ Q − l − mkn , (40)˜ Q lm,sgkn = i − ˜ Q lmkn + ˜ Q l − mkn − ˜ Q − lmkn + ˜ Q − l − mkn , (41)˜ Q lm,shkn = ˜ Q lmkn − ˜ Q l − mkn − ˜ Q − lmkn + ˜ Q − l − mkn , (42)for l = 0 , m >
0: ˜ Q m,qgkn = ˜ Q mkn + ˜ Q − mkn , (43)˜ Q m,sgkn = i − ˜ Q mkn + ˜ Q − mkn , (44)for l > , m = 0: ˜ Q l ,qgkn = ˜ Q l kn + ˜ Q − l kn , (45)˜ Q l ,qhkn = i( ˜ Q l kn − ˜ Q − l kn ) , (46)and for l = 0 , m = 0: ˜ Q ,qgkn = ˜ Q kn . (47)Note that both internal potentials, eqs. (31)-(32) and eqs. (36)-(37), depend onlyon the pre-calculated Q and inducing coefficients. Additionally, Q does not depend onlocation, making it particularly well-suited for satellite data. The methods presented in the previous sections enable estimation of time-series ofinducing coefficients in discrete non-overlapping time intervals (time windows). Let usdefine time intervals of length ∆ t . We assume that inducing coefficients are piece-wiseconstant within these time intervals. Then, convolution integrals such as eq. (20) or eqs. –8–anuscript submitted to JGR: Space Physics (28)-(29) can be approximated by discrete sums. For instance, for a time window cen-tered at t we can rewrite eq. (28) as g mn ( t ; σ ) ≈ N t (cid:88) j =0 I Q n ( j ; σ ) q mn ( t − j ∆ t ) , (48)where I Q n ( j ; σ ) = ˆ j ∆ t +∆ t/ j ∆ t − ∆ t/ Q n ( t ; σ )d t. (49)Similar expressions are obtained for other convolution integrals.With this, coefficients for a time window centered at t can be estimated by solv-ing a minimization problem, q ∗ , s ∗ = arg min q , s (cid:88) i ∈D t (cid:88) α ∈{ θ,φ } B oα,i − (cid:88) ( n,m ) ∈M + B mn,α ( (cid:126)r i , t ; σ ) , (50)where D t is a set of magnetic field observations in the current time window with B oα,i being the measured horizontal magnetic field component at location (cid:126)r i and time t i ; q , s ∈M + are vectors of inducing SH coefficients for the given time window; and the modelledfields are given by B mn,α ( (cid:126)r i , t ; σ ) = N t (cid:88) j =0 (cid:104) I m ( c ) n,α ( (cid:126)r i , j ; σ ) q mn ( t − j ∆ t ) + I m ( s ) n,α ( (cid:126)r i , j ; σ ) s mn ( t − j ∆ t ) (cid:105) . (51)For ground observations (see Section 2.4.1), we used I m ( c ) n,α ( (cid:126)r i , j ; σ ) = ˆ t j +∆ t/ t j − ∆ t/ B m ( c ) n,α ( (cid:126)r i , τ ; σ )d τ, (52) I m ( s ) n,α ( (cid:126)r i , j ; σ ) = ˆ t j +∆ t/ t j − ∆ t/ B m ( s ) n,α ( (cid:126)r i , τ ; σ )d τ. (53)Note that the two equations above are valid for magnetic fields computed either in 1-D or 3-D conductivity models.For satellite measurements (see Section 2.4.2) and a 1-D subsurface conductivitydistribution, we take I m ( c ) n,θ ( (cid:126)r i , j ; σ ) = − I Q n ( j ; σ ) (cid:18) ar i (cid:19) n +2 d P mn (cos θ )d θ (cid:12)(cid:12)(cid:12) θ = θ i cos( mφ i ) , (54) I m ( c ) n,φ ( (cid:126)r i , j ; σ ) = I Q n ( j ; σ ) (cid:18) ar i (cid:19) n +2 m sin θ i P mn (cos θ i ) sin( mφ i ) , (55) I m ( s ) n,θ ( (cid:126)r i , j ; σ ) = − I Q n ( j ; σ ) (cid:18) ar i (cid:19) n +2 d P mn (cos θ )d θ (cid:12)(cid:12)(cid:12) θ = θ i sin( mφ i ) , (56) I m ( s ) n,φ ( (cid:126)r i , j ; σ ) = − I Q n ( j ; σ ) (cid:18) ar i (cid:19) n +2 m sin θ i P mn (cos θ i ) cos( mφ i ) . (57)Similar, although more lengthy, expressions can be derived using eqs. (36)-(37) for a 3-D subsurface conductivity distribution.Note that since we have eliminated internal coefficients, it suffices to use only hor-izontal magnetic field components in eq. (50) to determine inducing coefficients. Thisallows for more accurate description of the inducing source since horizontal components –9–anuscript submitted to JGR: Space Physics are less sensitive to the currents induced in the subsurface compared to the vertical com-ponent (Kuvshinov, 2008). Since the problem is linear with respect to inducing coeffi-cients, we used a Huber-weighted robust regression method to find the minimizer of (50).For every time window, the performance of the model can be evaluated by meansof a R statistics, called coefficient of determination. To define it, let us assume that fora given field component all observations and modeled fields in a time window j are col-lected into vectors b obs j and b mod j such that r j = b obs j − b mod j (58)is the vector of residuals. Then R j = 1 − (cid:104) r j , r j (cid:105)(cid:104) b obs j − b obs j , b obs j − b obs j (cid:105) (59)is the coefficient of determination for time window j . Here b obs j denotes the mean valueof b obs j and (cid:104)· , ·(cid:105) is an inner product. Note that we assumed a uniform measurement er-ror of 1 nT when calculating R . Previous sections concentrated on evaluation of inducing coefficients. Once theyare estimated, we can evaluate induced coefficients that describe EM fields induced inthe planetary interior. This is useful in induction studies, where pairs of inducing andinduced coefficients are used to estimate subsurface transfer functions, which can be ul-timately inverted for the electrical conductivity distribution in the subsurface (P¨uthe &Kuvshinov, 2014).By adopting our approach, induced coefficients can be estimated from the radialcomponent alone. This is advantageous since the radial field exhibits higher sensitivityto the subsurface induction effects and was excluded from the estimation of the induc-ing coefficients (see eq. 50).Provided that the inducing coefficients q mn , s mn were estimated following the approachpresented in the previous section, the induced part of the total magnetic field can be iso-lated. In particular, for the observed radial magnetic field, B int , o r ( (cid:126)r, t ) = B o r ( (cid:126)r, t ) − B ext , o r ( (cid:126)r, t ) , (60)where B ext , o r ( (cid:126)r, t ) = − (cid:88) ( n,m ) ∈M + [ q mn ( t ) cos( mφ ) + s mn ( t ) sin( mφ )] n (cid:16) ra (cid:17) n − P mn (cos θ ) , (61)is the inducing part of the radial field (see eq. 23).Following eq. (25), the remaining induced part of the radial field above the groundcan be expanded as B int r ( (cid:126)r, t ) = (cid:88) k,l (cid:2) g lk ( t ) cos( kφ ) + h lk ( t ) sin( kφ ) (cid:3) ( k + 1) (cid:16) ar (cid:17) k +2 P lk (cos θ ) , (62)which is suitable for the estimation of the induced coefficients in a statistical manner.Specifically, we can estimate coefficients for a time bin centered at t = j ∆ t by solvinga minimization problem g ∗ , h ∗ = arg min g , h (cid:88) i ∈D t B int , o r,i − (cid:88) k,l (cid:2) g lk,j cos( kφ i ) + h lk,j sin( kφ i ) (cid:3) ( k + 1) (cid:18) ar i (cid:19) k +2 P lk (cos θ i ) . (63) –10–anuscript submitted to JGR: Space Physics
Therefore, by virtue of eqs. (50) and (63) pairs of inducing q ∗ , s ∗ and induced g ∗ , h ∗ coefficients can be estimated in time bins of constant length ∆ t , providing input datafor mantle conductivity studies. Note that estimation of inducing and induced coefficientscan be performed repeatedly with updated mantle conductivity models. We apply the developed methods to the ground geomagnetic observatory data. Specif-ically, we took a set of quality-controlled measurements of the hourly mean vector mag-netic field compiled by the British Geological Survey (Macmillan & Olsen, 2013). Weconcentrate here on the
Swarm era measurements by using data collected between 2013-12-01 and 2019-11-01. The model of the core and crustal fields as given by the Compre-hensive Inversion (CI) model (Sabaka et al., 2018) was subtracted. The distribution ofthe observatories over the time range used in this study is shown in Figure 1. We fur-ther excluded observatories poleward of the 56 ◦ and equatorward of 5 ◦ geomagnetic lat-itudes. Thus, the variations in the remaining data set are predominantly driven by themid latitude ionospheric and magnetospheric currents. The polar and equatorial latitudesare excluded because the present distribution of geomagnetic observatories can not ad-equately resolve spatiotemporal structures of the dominant current systems at these lat-itudes. We used nearly six years (2013-12-01 – 2019-11-01) of the geomagnetic field mea-surements taken by the
Swarm
Alpha and Bravo satellites. Similar to the observatorydata, core and crustal fields as given by the CI model were subtracted. The time win-dows of three hours were used, which corresponds to two full orbits and aims to improvethe data coverage within a window. Here, we concentrate on studying the EM inductioneffects of the large-scale magnetosphere currents of external origin and thus the day sidedata, namely between 5 AM and 7 PM local time, was excluded.
All transfer functions and corresponding impulse responses referred to as ”1-D” werecalculated by taking a conductivity model that consists of the 1-D conductivity profilefrom Grayver et al. (2017) with a 7000 S conductance layer that represents average con-ductance of the oceans and sediments. For the results referred to as ”3-D”, a laterallyheterogeneous conductivity shell of 1 / ◦ resolution was used to account for the varia-tions in the ocean bathymetry and thickness of sediments. For the 1-D case, transfer func-tion were calculated analytically, whereas 3-D transfer functions were calculated numer-ically by solving Maxwell’s equations in a spherical shell with a Finite Element code GoFEM(Grayver & Kolev, 2015; Grayver et al., 2019; Arndt et al., 2020).Figure 2 shows 1-D transfer functions and corresponding discrete impulse responses.As expected, we see that the decay rate for responses with higher degrees n is faster, im-plying that attenuation rate of the induced currents increases with the SH degrees of theinducing field. At periods of 1 year and longer, real part of the transfer function flattensas a result of the transient induction effect of the core, which has a finite conductivity(Vel´ımsk`y et al., 2003).Figure 3 shows a set of discrete impulse responses from the 3-D Q -matrix for dif-ferent external and internal degrees and orders. First of all, note that in 3-D, the ma- –11–anuscript submitted to JGR: Space Physics
Figure 1.
Top: Distribution of geomagnetic observatories. Location of geomagnetic observato-ries are denoted with circles. Filled circles show observatories used in this study, after discardinglocations at high and equatorial geomagnetic dipole latitudes. Bottom: number of used observa-tories over the time period of the study.
Figure 2.
Real (A) and imaginary (B) parts of the ˜ Q n transfer functions (eq. 27) for dif-ferent degrees n and 1-D conductivity profile of Grayver et al. (2017). The magnitudes of thecorresponding discrete impulse responses (eq. 49) are shown in plot (C).–12–anuscript submitted to JGR: Space Physics
Figure 3.
A selection of the 3-D discrete impulse responses from the Q lm,qgkn and Q lm,qhkn ma-trices (eqs. 36-37) due to the q (top row) and s (bottom row) inducing terms. Dashed linesdenote responses which are non zero only in the case of a 3-D conductivity distribution.–13–anuscript submitted to JGR: Space Physics
Figure 4.
Magnetic field discrete impulse responses (eq. 21) due to a q inducing field forthree magnetic field components (columns) at three locations: F¨urstenfeldbruck (FUR), Her-manus (HER) and Gan, Maldives (GAN). Both 1-D (dashed lines) and 3-D responses (solid lines)are shown. trix is dense, i.e. each inducing coefficients leads to infinitely many induced coefficients.However, we observe that the diagonal elements dominate the matrix, whereas off-diagonalentries are generally smaller and decay with the SH degree.Finally, Figure 4 shows examples for local impulse responses at several observatorylocations where both 1-D and 3-D responses are plotted to highlight the effect of the oceanand sedimentary cover on impulse responses. We see that the difference between 1-D and3-D responses is particularly large for island and coastal locations. We determined SH coefficients up to degree n max = 3 and order m max = 3 withinhourly time bins. The length of impulse responses was set to six months, thus transienteffects older than six months are neglected. This choice is justified since impulse responsesfor time lags larger than six months are ≤ − (see Figures 2-4), thus the transient ef-fects become negligible for the majority of practical applications. Other details pertainedto data pre-processing and the method of evaluating SH coefficients are given in Sections3.1 and 2.5, respectively.The coefficients were determined using both 1-D and 3-D impulse responses fromhorizontal magnetic field components ( B θ , B φ ). Subsequently, coefficient of determina-tion R was calculated for every time bin using eq. (59) and three components separately, –14–anuscript submitted to JGR: Space Physics
Figure 5.
Histograms of the R statistics (coefficient of determination) for individual mag-netic field components (from left to right: B r , B θ , B φ ) and all time windows. The R statisticswas determined following eq. (59) between observatory data and model predictions. The modeldetails are described in Section 4.2. including B r component, which was not used for the model construction. Figure 5 showshistograms of R coefficient for 1-D and 3-D models. One apparent observation is thesignificantly better fit of the radial component with a 3-D conductivity model. The fitfor horizontal components is virtually identical, the differences are minute and likely fallwithin the modelling and observation errors. Noteworthy that among all components,the highest coherency is observed for the longitudinal component. These observationsconfirm that our model, especially the one based on a 3-D conductivity model, has a pre-dictive power.To test how the model performs during different magnetic conditions, we furtherplot histograms of R statistics for times when magnetic variations are dominated bymagnetospheric disturbances (here defined as | Dst | >
40 nT) in Figure 6. Althoughwe still observe a significant improvement in coherency for B r , generally the correlationis lower for B r , B θ , whereas it remains high for the longitudinal component. Further, weplot R histograms for times when Kp ≤
2. The reason to use Kp instead of Dst thistime is to emphasize the quiet ionosphere conditions. Similar to the examples with dis-turbed magnetosphere, we again observe significant improvements in the radial compo-nent for the 3-D model. In comparison with the previous case, however, we see system-atically higher R values for all components. Therefore, our model exhibits a better fitduring quiet times. Further, the improved fit of the B r enabled by the 3-D model is ob-served for all times and magnetic conditions, indicating that proper inclusion of the oceaneffect is essential when modeling both magnetospheric and ionospheric variations.To better quantify the effect of the improved fit due to the usage of a 3-D model,we calculated the ratio of 3-D and 1-D R values for radial magnetic field componentat all observatory locations. These values, plotted as a function of the distance to theshoreline, are shown in Figure 8. We observe the improved fit at virtually all locationswith the most significant improvement up to a factor of 11 for observations that are ≤
200 km from the coast. However, even locations as far as 3000 km exhibit considerablybetter fit. This is explained by including the conductance of continental sediments in our3-D model.Finally, we inspect the observed and modeled time series at a selection of coastaland island observatories. Here, we also added predictions based on the Dst index, cal- –15–anuscript submitted to
JGR: Space Physics
Figure 6.
Same as Figure 5, but restricted to time windows when | Dst | >
40 nT.
Figure 7.
Same as Figure 5, but restricted to time windows when Kp ≤ Figure 8.
Ratio of the 3-D to 1-D models R coefficients for B r field at individual observato-ries plotted versus distance to the shoreline. Values larger than one indicate improvement overthe 1-D model. –16–anuscript submitted to JGR: Space Physics culated as B Dst r ( (cid:126)r, t ) = (cid:2) Est( t ) − t ) (cid:3) cos θ (64) B Dst θ ( (cid:126)r, t ) = − (cid:2) Est( t ) + Ist( t ) (cid:3) sin θ (65) B Dst φ ( (cid:126)r, t ) = 0 , (66)where Dst( t ) = Est( t ) + Ist( t ) is a sum of inducing and induced terms (Maus & Wei-delt, 2004; Olsen et al., 2005). Figures 9-10 each show one week of observed magneticfield variations and model predictions. These periods were chosen since they cover bothmagnetically disturbed and quiet conditions.The origin of the discrepancy in amplitude between the observed and modelled fieldsis twofold: (i) we used the global average mantle conductivity profile whereas in realitythe bulk subsurface conductivity varies laterally; and (ii) slightly larger discrepancy forquiet times indicates that the adopted SH parameterization with n max = 3 , m max =3 is still insufficient to explain these variations, mostly related to ionospheric currents(Guzavina et al., 2019; Schmucker, 1999).Additionally, Figures 11-12 show spatial distribution of the magnetic field as pre-dicted by estimated external coefficients and a 3-D conductivity model. Much strongerinfluence of 3-D EM induction effects in the B r components are clearly visible. Most ofthese effects occur near coastal areas and strong lateral conductivity gradients. Swarm ob-servations
In this section, the model of inducing coefficients was determined by using only satel-lite data, which was described in Section 3.2. Since we work with night-side data andtwo satellites, we determined SH coefficients up to degree n max = 2 and order m max =1 using time bins of 3 hours. Therefore, the resolution of this model is much lower thanthe model in previous section that was based on observatory data. Other parameters per-tained to data pre-processing and evaluation of SH coefficients are described in Sections3.2 and 2.5, respectively.As in the previous section, we first look at the distribution of R statistics for alltime bins and magnetic field components (see Figure 13). First observation that we makeis that R values are very similar between 1-D and 3-D models, indicating that 3-D in-duction effect from the ocean is largely attenuated at satellite altitudes. Interestinglythat now we also have much higher R values for the radial component compared to the B θ , even though B r was not used in the construction of the model. One possible expla-nation are signals that mostly affect horizontal ( B θ , B φ ) components at mid latitudes,such as those generated by F-region ionospheric currents (Olsen, 1997). These signalscannot be explained by our low-resolution parameterization that is based on the poten-tial field assumption. To test this hypothesis, histograms limited to the time windowsfor which Kp ≤ R for the B θ component.Finally, Figure 15 plots time series of the q coefficient determined using the ob-servatory and satellite data. For reference, we also plot the D st index. We observe verygood match between coefficients estimated from satellite and observatory data, confirm-ing the validity of both models and approaches. The EM induction effect from a time-varying magnetic field significantly influencesmagnetic field observations, where it can be both a polluting signal to be removed or aprimary signal to study (e.g. mantle induction and space weather applications). We showed –17–anuscript submitted to
JGR: Space Physics
Figure 9.
Time series of observed and modelled variations in horizontal ( δX = − δB θ ) andradial ( δZ = − δB r ) components at a set of observatories, ordered by latitude. Predictions basedon 1-D and 3-D conductivity models are shown along with D st -based fields (eq. 64). The offsetbetween dotted lines is 100 nT. Lower panels show corresponding D st and Kp indices. See Figure1 for locations of selected observatories. –18–anuscript submitted to JGR: Space Physics
Figure 10.
Same as Figure 9, but for a different time period.–19–anuscript submitted to
JGR: Space Physics
Figure 11.
Maps of the radial (top) and horizontal (bottom) components of modelled mag-netic field variations at a surface. Predictions based on the 3-D conductivity model for a givenUT instance are shown. The recorded D st index value at this instance was -155 nT. Significantcoastal EM induction effects are visible in the radial component.–20–anuscript submitted to JGR: Space Physics
Figure 12.
Same as Figure 11, but for a different time. The D st index value at this instancewas -95 nT. –21–anuscript submitted to JGR: Space Physics
Figure 13.
Histograms of the R statistics (coefficient of determination) for individual mag-netic field components and all time windows. The R statistics was determined following eq. (59)between observatory data and predictions from satellite data model described in Section 4.3. Figure 14.
Same as Figure 13, but restricted to time windows when Kp ≤ JGR: Space Physics
Figure 15.
Time series of the first zonal SH coefficient, q , as given by satellite (Section 4.3)and observatory (Section 4.2) data based models in geomagnetic coordinate frame. Two fivemonths intervals featuring quiet and disturbed magnetic conditions are shown. For comparison,the negative D st magnetic index is plotted. Systematic offset in D st against q seen in Figures 15is due to the absence of stable quiet time ring current in the D st index.–23–anuscript submitted to JGR: Space Physics that the inducing currents of ionospheric and magnetospheric origin can be effectivelyestimated while the effect of the planetary induced response is modeled. This work haspresented a unified framework for modeling EM induction effects in ground and satel-lite data by means of time domain impulse responses due to arbitrary external sourcesand in presence of a 3-D subsurface conductivity distribution. This approach is amenableto integrate with models that involve constantly augmented time series and require ”onthe fly” updates of geomagnetic models.We have elaborated the underlying mathematical machinery for the case when ba-sis functions used for spatial parameterization of magnetic field are given by sphericalharmonic functions. This choice was made owing to the ubiquity of SH basis in Earth’sand planetary magnetism community. However, the approach is general and straight-forward to extend to other basis functions should practical applications demand this.We further showed that the effects from heterogeneity in subsurface electrical con-ductivity can dominate the radial magnetic field component and should be accountedfor provided that some knowledge about 3-D subsurface conductivity structure for Earthis available. Contrary to the common presumption, the 3-D effects are significant dur-ing both quiet and disturbed magnetic conditions since the induction effect is transient,hence widely used selection criteria based on instant values of magnetic indices and lo-cal time can not completely eliminate the effects of EM induction, rendering the mod-eling approach presented here a suitable alternative that accounts for its transient na-ture.
Appendix A Properties of transfer functions and impulse responses
Convolution integrals such as in eqs. (20) and (28)-(29) represent a response of amedium to a time-varying extraneous current. These relations follow from the (often omit-ted) properties of a physical system that we model. We state these properties here anddiscuss implications. Our presentation closely follows a more detailed analysis by Svetov(1991).1.
Linearity allows us to define a response, ζ ( t ), of a medium at time t to an ex-traneous forcing as ζ ( t ) = ˆ ∞−∞ F ( t, t (cid:48) ) χ ( t (cid:48) ) dt (cid:48) , (A1)where χ is the extraneous forcing that depends on time t (cid:48) and F ( t, t (cid:48) ) is the mediumGreen’s function that does not depend on the amplitude of the exerted force.2. Stationarity implies that the response of a medium does not depend on the timeof occurrence of the excitation. In this case F ( t, t (cid:48) ) ≡ f ( t − t (cid:48) ) and eq. (A1) canbe rewritten as a convolution integral ζ ( t ) = ˆ ∞−∞ f ( t − τ ) χ ( τ ) dτ = ˆ ∞−∞ f ( τ ) χ ( t − τ ) dτ, (A2)where f ( t ) represents the impulse response of a medium. In frequency domain,the convolution integral reduces to˜ ζ ( ω ) = ˜ f ( ω ) ˜ χ ( ω ) , (A3)where ˜ f ( ω ) is called the transfer function and we used tilde sign (˜ · ) to denote com-plex quantities. Eqs. (A2) and (A3) are related through the Fourier transform˜ f ( ω ) = ˆ ∞−∞ f ( t ) e i ωt d t. (A4)3. Since we work in time domain with a real valued forcing, the impulse response isalso real . To see implications of this, let us define the inverse Fourier transform –24–anuscript submitted to JGR: Space Physics of ˜ f ( ω ) = f R ( ω ) + i f I ( ω ) as f ( t ) = 12 π ˆ ∞−∞ ˜ f ( ω ) e − i ωt d ω = 12 π ˆ ∞−∞ [ f R ( ω ) cos( ωt ) + f I ( ω ) sin( ωt )] d ω + i2 π ˆ ∞−∞ [ f I ( ω ) cos( ωt ) − f R ( ω ) sin( ωt )] d ω. (A5)For an impulse response to be real, the last term in the integral (A5) has to van-ish. This is possible only if f R ( ω ) and f I ( ω ) are even and odd functions of the an-gular frequency ω , respectively. Therefore, eq. (A5) reduces to f ( t ) = 1 π ˆ ∞ [ f R ( ω ) cos( ωt ) + f I ( ω ) sin( ωt )] d ω. (A6)4. Impulse response is causal . This property implies that f ( t ) = 0 for t <
0. Un-der this assumption, the convolution integral (A2) can be recast to ζ ( t ) = ˆ ∞ f ( τ ) χ ( t − τ ) dτ = ˆ t −∞ f ( t − τ ) χ ( τ ) dτ. (A7)Due to causality and taking into account eq. (A6), the impulse response can bedetermined by using either only real or imaginary part of ˜ f ( ω ): f ( t ) = 2 π ˆ ∞ f R ( ω ) cos ( ωt )d ω (A8)= − π ˆ ∞ f I ( ω ) sin ( ωt )d ω. (A9)Note that for the sake of clarity, dependence on spatial variables and electrical con-ductivity pertinent to our application was omitted from the equations above.In practice, we observed that using sine transform (A9) results in a slightly bet-ter accuracy compared to the cosine transform (A8) given the same filter length. Appendix B Digital Linear Filters
In order to carry out the sine transform (A9) efficiently, we applied the linear dig-ital filter (DLF) method. DLF was introduced to geophysics by Ghosh in the early 70s(Ghosh, 1971a, 1971b), as a means of fast computations for geoelectric resistivity responses.The method was subsequently improved and expanded to other methods by many au-thors, and a lot of filters have been published. There are two particular developments,out of all these improvements, which are relevant for our application: (1) The kernel un-der consideration were early on always Bessel functions of some sort, and it was Anderson(1973) who first applied it to Fourier sine and cosine transforms. (2) If the kernel com-putation is very expensive the lagged-convolution type DLF introduced by Anderson (1975)is very powerful, as additional times come at no or very little extra cost due to the reuseof the already computed kernels for new times. Although the use of DLF in geophysicsis focused on Hankel and Fourier transforms in electromagnetics, the method itself worksfor any linear transform.Werthm¨uller et al. (2019) presented a tool to design filters for any linear transformprovided that there exist (a) an analytical transform pair, or (b) a numerical computa-tion in both domains with sufficient accuracy and precision over a wide range of argu-ment values. We refer to that publication for an in-depth review of DLF in geophysics.Using substitutions t = e x and ω = e − y and multiplying by e x we can rewrite(A9) as a convolution integral and approximate it by a N -point digital filter η as (Anderson, –25–anuscript submitted to JGR: Space Physics f ( t ) ≈ − π N (cid:88) n =1 f I ( b n /t ) η n t , (B1)where the log-spaced filter abscissa values b n are a function of spacing ∆ and shift ν , b n (∆ , ν ) = exp [∆( −(cid:98) ( N + 1) / (cid:99) + n ) + ν ] . (B2)The optimal values for η n , ∆ and ν in eqs. (B1)-(B2) were found by following themethod of Werthm¨uller et al. (2019). In this work, we designed a 50-point filter such thatit requires as few values of ˜ f ( ω ) as possible without compromising accuracy. To this end,we used the following analytic transform pair π exp ( − ab )2 = ˆ ∞ xa + x sin( xb ) dx. (B3)The Figure B1 shows the designed filter and its performance for the chosen analytic pair. Figure B1.
Left: minimum recovered value of the analytic pair (B3) as a function of spac-ing and shift. Center: filter values for the best filter with ∆ = 0 .
114 and ν = 1 .
07. Right: theperformance of the filter on the eq. (B3).
Note that the naive application of eq. (B1) will require calculating ˜ f ( ω ) at N × N t frequencies, where N t is the length of an impulse response in time domain. This num-ber can be drastically reduced by invoking the aforementioned lagged convolution ap-proach. To give an example, for the one year long impulse response with the hourly timestep (i.e., N t = 8766) our filter required evaluating a maximum of 112 frequencies thatrange ≈
12 decades.
Acknowledgments
This work was supported by the ESA through the Swarm DISC project. The staff of thegeomagnetic observatories and INTERMAGNET are thanked for supplying high-qualityobservatory data, and BGS are thanked for making quality-controlled observatory hourlymean values openly available. Satellite and observatory input data used in this studyas well as estimated time series of SH coefficients can be retrieved from https://doi.org/10.5281/zenodo.4047833 . Subsurface conductivity model is available at https://doi.org/10.5281/zenodo.4058852 . We thank Chris Finlay, Nils Olsen and Jakub Vel´ımsk´yfor insightful discussions that facilitated this study.
References
Anderson, W. L. (1973).
Fortran IV programs for the determination of the transienttangential electric field and vertical magnetic field about a vertical magnetic –26–anuscript submitted to
JGR: Space Physics
Table B1.
Filter used in this study to approximate the sine transform A9. base value0.020351539057584585 0.0156125158035318530.02394618633584547 -0.094629034117499140.028175748203049585 0.31177704934884790.0331523682171173 -0.73317330329696880.039007997604279074 1.38871956070486950.045897893843668 -2.2566570665789280.05400473719916161 3.28053589354546160.06354347434512761 -4.37800141257025550.0747670175110597 5.4726339812070250.08797295025351011 -6.4876800553071040.10351141765367007 7.3792030815479220.121794410143077 -8.1014066993795840.1433066871108987 8.6557110191325680.16861862992378368 -9.0176101064279290.19840136514614606 9.2294905854077830.23344455894136176 -9.2632114693298430.2746773544586444 9.2016250182173830.3231930073442032 -8.9810174201598620.38027787256818374 8.7411935867071280.44744551113066483 -8.3396369111785980.5264768209596123 8.0172526943391630.6194672560404727 -7.4843852587480440.728882385756068 7.1689394846408280.8576232675496689 -6.52052834371884151.009103366216787 6.2983821480117871.1873390592811512 -5.5178357294908351.3970561281348315 5.4673762615898121.6438150584726328 -4.536620002704381.9341584722647591 4.67312136269356952.2757846003123405 -3.6912900547884412.6777513948763136 3.7823710985624273.1507166942679676 -3.22814399338347753.7072208071793002 2.55925464376196574.36201900925792 -3.2321196658135575.1324727678156865 1.51263466215251446.039010067691346 -2.24901577981362837.105667043450954 2.6249173695595768.360725278884502 -0.387783447640519759.837461671301405 0.769756567858555611.575030742696057 -2.860710788307397213.619502791580832 3.119988106045642716.025085411278443 -1.912035359909057418.855560762285485 0.80975332601978922.185976706870886 -0.264192336906859126.104636645033704 0.0708812213054222130.715440810780017 -0.01612140316180455536.14064110637662 0.003087124452875387542.524082516878885 -0.0004729822575874928350.03501704852909 5.153158194873856e-0558.87259131464845 -2.9723596837103312e-06 –27–anuscript submitted to
JGR: Space Physics dipole for an m-layer stratified earth by numerical integration and digital linearfiltering (Vol. PB221240; Tech. Rep.). US Geological Survey.Anderson, W. L. (1975).
Improved digital filters for evaluating fourier and hankeltransform integrals (Tech. Rep.). US Geological Survey.Arndt, D., Bangerth, W., Blais, B., Clevenger, T. C., Fehling, M., Grayver,A. V., . . . Wells, D. (2020). The deal.ii library, version 9.2.
Jour-nal of Numerical Mathematics (0), 000010151520200043. Retrieved from doi: https://doi.org/10.1515/jnma-2020-0043Chulliat, A., Vigneron, P., & Hulot, G. (2016). First results from the Swarm dedi-cated ionospheric field inversion chain.
Earth, Planets and Space , (1), 1–18.Gauss, C. F. (1877). Allgemeine Theorie des Erdmagnetismus. In Werke (pp. 119–193). Springer.Ghosh, D. P. (1970).
The application of linear filter theory to the direct interpre-tation of geoelectrical resistivity measurements (Ph.D. Thesis, TU Delft). doi:http://resolver.tudelft.nl/uuid:88a568bb-ebee-4d7b-92df-6639b42da2b2Ghosh, D. P. (1971a). The application of linear filter theory to the direct interpre-tation of geoelectrical resistivity sounding measurements.
Geophysical Prospect-ing , (2), 192–217. doi: 10.1111/j.1365-2478.1971.tb00593.xGhosh, D. P. (1971b). Inverse filter coefficients for the computation of apparentresistivity standard curves for a horizontally stratified earth. GeophysicalProspecting , (4), 769–775. doi: doi:10.1111/j.1365-2478.1971.tb00915.xGrayver, A. V., & Kolev, T. V. (2015). Large-scale 3D geoelectromagnetic modelingusing parallel adaptive high-order finite element method. Geophysics , (6),E277–E291.Grayver, A. V., Munch, F. D., Kuvshinov, A. V., Khan, A., Sabaka, T. J., &Tøffner-Clausen, L. (2017). Joint inversion of satellite-detected tidal andmagnetospheric signals constrains electrical conductivity and water content ofthe upper mantle and transition zone. Geophysical research letters , (12),6074–6081.Grayver, A. V., van Driel, M., & Kuvshinov, A. V. (2019). Three-dimensional mag-netotelluric modelling in spherical earth. Geophysical Journal International , (1), 532–557.Guzavina, M., Grayver, A., & Kuvshinov, A. (2019). Probing upper mantle electri-cal conductivity with daily magnetic variations using global-to-local transferfunctions. Geophysical Journal International , (3), 2125–2147.Honkonen, I., Kuvshinov, A., Rast¨atter, L., & Pulkkinen, A. (2018). Predictingglobal ground geoelectric field with coupled geospace and three-dimensionalgeomagnetic induction models. Space Weather , (8), 1028–1041.Kelbert, A. (2020). The role of global/regional earth conductivity models in naturalgeomagnetic hazard mitigation. Surveys in Geophysics , (1), 115–166.Koch, S., & Kuvshinov, A. (2013). Global 3-D EM inversion of Sq variations basedon simultaneous source and conductivity determination: concept validationand resolution studies. Geophysical Journal International , (1), 98–116.Kuvshinov, A. (2008). 3-D global induction in the oceans and solid Earth: Recentprogress in modeling magnetic and electric fields from sources of magneto-spheric, ionospheric and oceanic origin. Surv Geophys , , 139-186. doi:10.1007/s10712-008-9045-zKuvshinov, A., Avdeev, D., & Pankratov, O. (1999). Global induction by sq anddst sources in the presence of oceans: bimodal solutions for non-uniform spher-ical surface shells above radially symmetric earth models in comparison toobservations. Geophysical journal international , (3), 630–650.Macmillan, S., & Olsen, N. (2013). Observatory data and the swarm mission. Earth,Planets and Space , (11), 15. –28–anuscript submitted to JGR: Space Physics
Maus, S., & Weidelt, P. (2004). Separating the magnetospheric disturbance magneticfield into external and transient internal contributions using a 1d conductivitymodel of the earth.
Geophysical research letters , (12).Munch, F. D., Grayver, A. V., Guzavina, M., Kuvshinov, A. V., & Khan, A. (2020).Joint inversion of daily and long-period geomagnetic transfer functions re-veals lateral variations in mantle water content. Geophysical Research Letters ,e2020GL087222.Olsen, N. (1997). Ionospheric f region currents at middle and low latitudes estimatedfrom magsat data.
Journal of Geophysical Research: Space Physics , (A3),4563–4576.Olsen, N. (1999). Induction studies with satellite data. Surveys in Geophysics , (3-4), 309–340.Olsen, N., Glassmeier, K.-H., & Jia, X. (2010). Separation of the magnetic field intoexternal and internal parts. Space science reviews , (1-4), 135–157.Olsen, N., & Kuvshinov, A. (2004). Modeling the ocean effect of geomagneticstorms. Earth, planets and space , (5), 525–530.Olsen, N., Sabaka, T. J., & Lowes, F. (2005). New parameterization of external andinduced fields in geomagnetic field modeling, and a candidate model for IGRF2005. Earth, planets and space , (12), 1141–1149.Price, A. (1967). Electromagnetic induction within the Earth. In International geo-physics (Vol. 11, pp. 235–298). Elsevier.P¨uthe, C., & Kuvshinov, A. (2014). Mapping 3-d mantle electrical conductivity fromspace: a new 3-d inversion scheme based on analysis of matrix q-responses.
Geophysical Journal International , , 768-784. doi: 10.1093/gji/ggu027P¨uthe, C., Manoj, C., & Kuvshinov, A. (2014). Reproducing electric field ob-servations during magnetic storms by means of rigorous 3-d modelling anddistortion matrix co-estimation. Earth, Planets and Space , (1), 162.Sabaka, T. J., Olsen, N., & Purucker, M. E. (2004). Extending comprehensive mod-els of the Earth’s magnetic field with ørsted and CHAMP data. GeophysicalJournal International , (2), 521–547.Sabaka, T. J., Olsen, N., Tyler, R. H., & Kuvshinov, A. (2015). CM5, a pre-Swarmcomprehensive geomagnetic field model derived from over 12 yr of CHAMP,ørsted, SAC-C and observatory data. Geophysical Journal International , (3), 1596–1626.Sabaka, T. J., Tøffner-Clausen, L., Olsen, N., & Finlay, C. C. (2018). A compre-hensive model of Earth’s magnetic field determined from 4 years of Swarmsatellite observations. Earth, Planets and Space , (1), 130.Schmucker, U. (1985). Magnetic and electric fields due to electromagnetic induc-tion by external sources, electrical properties of the Earth’s interior. Landolt-Bornstein , New–Series.Schmucker, U. (1999). A spherical harmonic analysis of solar daily variations in theyears 1964–1965: response estimates and source fields for global induction—ii.results.
Geophysical Journal International , (2), 455–476.Schuster, A. (1889). XV. the diurnal variation of terrestrial magnetism. Philosophi-cal Transactions of the Royal Society of London.(A.) (180), 467–518.Sun, J., & Egbert, G. D. (2012). Spherical decomposition of electromagnetic fieldsgenerated by quasi-static currents.
GEM-International Journal on Geomathe-matics , (2), 279–295.Sun, J., Kelbert, A., & Egbert, G. D. (2015). Ionospheric current source modelingand global geomagnetic induction using ground geomagnetic observatory data. Journal of Geophysical Research: Solid Earth , (10), 6771–6796.Svetov, B. S. (1991). Transfer functions of the electromagnetic field (in russian). Fizika Zemli (1), 119–128.Thomson, A. W., & Lesur, V. (2007). An improved geomagnetic data selection al-gorithm for global geomagnetic field modelling.
Geophysical Journal Interna- –29–anuscript submitted to
JGR: Space Physics tional , (3), 951–963.Vel´ımsk´y, J., Everett, M. E., & Martinec, Z. (2003). The transient Dst electro-magnetic induction signal at satellite altitudes for a realistic 3-D electricalconductivity in the crust and mantle. Geophysical research letters , (7).Vel´ımsk´y, J., & Martinec, Z. (2005). Time-domain, spherical harmonic-finite elementapproach to transient three-dimensional geomagnetic induction in a sphericalheterogeneous Earth. Geophysical Journal International , (1), 81–101.Werthm¨uller, D., Key, K., & Slob, E. C. (2019). A tool for designing digital fil-ters for the hankel and fourier transforms in potential, diffusive, and wavefieldmodeling. Geophysics , (2), F47–F56.Yamazaki, Y., & Maute, A. (2017). Sq and EEJ—A review on the daily variationof the geomagnetic field caused by ionospheric dynamo currents. Space ScienceReviews , (1-4), 299–405.(1-4), 299–405.