Velimsky_PhD

Charles University in Prague Faculty of Mathematics and Physics ELECTROMAGNETIC INDUCTION IN A HETEROGENEOUS EARTH’S MANTLE: TIME-DOMAIN MODELLING Jakub Vel ´ ımsk´ y Prague 2003 Charles University in Prague Faculty of Mathematics and Physics ELECTROMAGNETIC INDUCTION IN A HETEROGENEOUS EARTH’S MANTLE: TIME-DOMAIN MODELLING Jakub Vel ´ ımsk´ y Ph.D. thesis Supervisor: Zdenˇ ek Martinec Department of Geophysics V Holeˇ soviˇ ck´ ach 2 180 00 Praha 8 Czech Republic m´ amˇ e a t´ atovi Dˇ ekuji sv´ emu ˇ skoliteli Zdeˇ nku Martincovi za veˇ skerou pomoc, radu, motivaci i kritiku, kter´ ych se mi od nˇ ej v pr˚ ubˇ ehu m´ eho studia hojnˇ e dost´ avalo. Dˇ ekuji tak´ e Jiˇ r ´ ımu Zahradn ´ ıkovi, Ctiradu Matyskovi, Ladislavu Hanykovi i ostatn ´ ım koleg˚ um z Katedry geofyziky MFF UK za b´ ajeˇ cnou atmosf´ eru, kter´ a zde panuje. I thank Mark Everett for his help, motivation, and for fruitful discussions we had during my staying at Texas A&M. Contents 1 Introduction 5 1.1 EM induction as a tool for Earth’s mantle studies . . . . . . . . . . . . . 5 1.2 Overview of previous 3-D approaches . . . . . . . . . . . . . . . . . . . . 6 2 Time-domain, spectral-finite element approach 9 2.1 Formulation of the initial-value problem . . . . . . . . . . . . . . . . . . 9 2.1.1 Classical formulation . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.1.2 Integral (weak) formulation . . . . . . . . . . . . . . . . . . . . . 14 2.1.3 Time integration scheme . . . . . . . . . . . . . . . . . . . . . . . 16 2.1.4 Spherical harmonic-finite element parameterization . . . . . . . . 18 2.2 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 2.2.1 Nested-sphere conductivity models . . . . . . . . . . . . . . . . . 25 2.2.2 Excitation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 2.2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3 Time-domain, 3-D finite element approach 35 3.1A−Φ,Uformulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.1.1 Classical formulation . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.1.2 Integral (weak) formulation . . . . . . . . . . . . . . . . . . . . . 40 3.1.3 Time integration scheme . . . . . . . . . . . . . . . . . . . . . . . 41 3.1.4 Discretization using 3-D nodal finite elements . . . . . . . . . . . 43 3.2 Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.2.1 Axially symmetric conductivity model . . . . . . . . . . . . . . . 47 3.2.2 Nested-sphere conductivity model . . . . . . . . . . . . . . . . . . 49 4 The transient Dst induction signal at satellite altitudes 55 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.2 Conductivity model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 4.3 Exciting field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 5 Conclusions 67 References 69 4 A Spherical harmonic functions 73 B Existence and uniqueness of the integral solution 77 C Tetrahedral mesh 81 D Nodal finite element functions 89 E Moving least squares interpolation 93 F Notation 95 F.1 Used symbols . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 F.2 Overview of functional spaces . . . . . . . . . . . . . . . . . . . . . . . . 96 5 Chapter 1 Introduction 1.1 EM induction as a tool for Earth’s mantle studies Studies of the electrical conductivity are an important constraint on the physical prop- erties of the Earth’s mantle. The electrical conductivity is particularly sensitive to tem- perature, chemical composition of mantle rocks and minerals, and to the presence of interstitial fluids in the form of aqueous phases or partial melts. Our knowledge of the conductivity distribution complements the information provided by seismology, simula- tions of thermochemical mantle convection, geoid studies, laboratory measurements, and other geophysical disciplines in the efforts to understand the physical processes in the mantle. The response of the conductive Earth’s mantle to long-periodic variations of the iono- spheric and magnetospheric currents represents the most effective way to determine the mantle conductivity and it has been intensively studied for a long time. The first inver- sions of the geomagnetic variations measured at the surface assumed a spherically sym- metric distribution of conductivity in the Earth’s mantle (Price, 1930; Lahiri & Price, 1939; Rikitake, 1950; Banks, 1969). Results obtained by mantle convection modelling and by the 3-D seismic tomography suggest that the lateral variations of temperature in the mantle are of orders of 102–103K on the scales of 102–103km (Duffy & Hemley, 1995). These results and the laboratory measurements of conductivity under mantle pressure and temperature conditions (e.g. Shankland et al., 1993; Xu et al., 1998, 2000) yield that the lateral conductivity varia- tions of orders of magnitude can be expected in the mantle. Local 1-D interpretations of carefully processed data strongly support this idea. Schultz & Larsen (1987) demon- strated that no single 1-D conductivity model is compatible with geomagnetic responses observed at geographically distributed observatories. The existence of lateral conduc- tivity variations in the upper and mid-mantle has been confirmed in other local 1-D and regional studies (Roberts, 1984; Schultz & Larsen, 1990; Schultz, 1990; Petersons & Anderssen, 1990; Tarits, 1994; Olsen, 1998; Semenov, 1998, a.o.). These results represented a strong motivation for the development of efficient forward solvers of the electromagnetic (EM) induction problem in a radially and laterally inho- mogeneous Earth’s mantle. A brief overview of these solvers using various mathematical approaches and numerical methods is given in the next section. 6 1. INTRODUCTION One of the main obstacles in the inversion of the ground based long-periodic data to reveal a 3-D conductivity structure is the uneven spatial distribution of the geomag- netic observatories, especially the poor coverage of the oceanic areas. On the other hand, the satellite magnetometer data such as those obtained from the MAGSAT, Ørsted and CHAMP spacecrafts have excellent spatial coverage and offer an unprecedented oppor- tunity to improve our knowledge of electrical conductivity in the Earth’s mantle. The problem of EM induction in the mantle is traditionally approached in the fre- quency domain. The equation of EM induction is solved with the assumption of harmonic signals at discrete frequencies usually from the range of 10−2–100cpd, covering the spec- trum that corresponds to solar quiet (Sq), and disturbed (Dst) geomagnetic variations. However, the transient signals, such as magnetic storms, are sources of strong induction effects in the mantle due to abrupt changes of the system of external inducing electrical currents in the ionosphere and magnetosphere. It is inconvenient to study the Earth’s response to transient signals in the frequency domain because they have a broad spec- trum. Moreover, the complicated spatio-temporal distribution of satellite data favors a time-domain approach. In this work we present two different time-domain approaches to the forward EM induction problem in a heterogeneous mantle. The spectral-finite element method intro- duced in Chapter 2 uses a traditional parameterization by spherical harmonic function. Chapter 3 introduces an approach based on the 3-D finite element discretization, a method common in the world of electrical engineering. In Chapter 4 which is the extended ver- sion of a GRL paper (Vel´ ımsk´ y et al., 2003) we employ the former method to compute the response of a realistic 3-D model of mantle conductivity to a transient signal with simplified spatial structure and realistic time evolution. 1.2 Overview of previous 3-D approaches The development of forward solvers of the EM induction problem in a spherical Earth with laterally varying conductivity has been subject of active research recently. Following the traditional frequency-domain approach, Fainberg et al. (1990) assumed a thin het- erogeneous conductive sheet representing the high contrast between conducting oceans and resistive continents with underlying radially stratified medium. The problem was formulated as an integral equation for surface currents in the sheet and solved by means of an iterative-dissipative method. The method was later generalized for an arbitrary position of the heterogeneous thin sheet in the spherically symmetric Earth (Kuvshinov & Pankratov, 1994). Zhang & Schultz (1992) developed a perturbation method for a radially symmetric conductivity model superposed by small lateral heterogeneities. The magnetic field was expressed in terms of toroidal and poloidal modes which were expanded into perturbation 1.2. OVERVIEW OF PREVIOUS 3-D APPROACHES 7 series truncated at 2. The zero-th, first, and second order terms were parameterized by spherical harmonic functions and finite differencing was used to solve the ordinary differential equations with respect to radius. Everett & Schultz (1996) applied the finite element method to the problem formulated in terms of electromagnetic potentials in the mantle and surrounding insulating atmo- sphere. The tetrahedralization of the sphere was based on Delauney triangulation of a spherical surface, piecewise linear nodal finite elements were used to parameterize the potentials, and iterative refinement scheme with incomplete LU factorization was used to solve the sparse linear system. Martinec (1999a) formulated the problem in the integral sense for the magnetic induc- tion vector B. Using vector spherical harmonic functions and piecewise linear elements for angular and radial parameterization, respectively, he constructed a complex sparse linear system that was solved by biconjugate gradient iteration. An efficient precondi- tioner was based on a spherically symmetric conductivity profile. The original formulation prescribed the boundary condition on the Earth’s surface in terms of total toroidal elec- trical intensity. However, the method was later modified in such a way that the Gauss coefficient of the external inducing field can be prescribed on the Earth’s surface. A staggered-grid finite difference method was implemented by Uyeshima & Schultz (2000). The computational domain was divided into curved rectangular H-prisms with a constant conductivity value assigned to each one. The integral form of the Amp` ere law for magnetic intensity Hwas cast in each H-prism, while the integral form of the Faraday law for electrical intensity Ewas cast in staggered E-prisms whose vertices were positioned in centers of the H-prisms. The sparse linear system was solved using the minimum residual conjugate gradient algorithm with incomplete Cholesky preconditioning. The computational domain enclosed the conductive mantle as well as the atmosphere which was approximated by prescribing a small conductivity value. An edge-based finite element formulation was introduced by Yoshimura & Oshiman (2002). This parameterization allowed to formulate the problem in terms of only one vector magnetic potential. The time-domain approach to the EM induction problem was recently applied by Hamano (2002). The problem was formulated in terms of poloidal and toroidal functions. Spherical harmonic expansion and finite differences were used in the angular and radial parameterization, respectively. The Crank-Nicolson scheme was employed for the time integration. Concurrently with the development of the 3-D methods, the question of their testing and validation was addressed. Everett & Schultz (1995) introduced a semi-analytical frequency-domain solution for a conductivity model consisting of a homogeneous host sphere and a homogeneous spherical inclusion nested eccentrically on the axis of sym- metryz. Martinec (1998) generalized the solution to the case of multiple off-axis sub- 8 1. INTRODUCTION sequently nested spheres. The spectral-finite element approach to the axially symmetric 2-D EM induction was implemented by Martinec (1997) and Martinec et al. (2002) in the frequency-domain and time-domain, respectively. 9 Chapter 2 Time-domain, spectral-finite element approach 2.1 Formulation of the initial-value problem 2.1.1 Classical formulation Let us start with notation and basic assumptions. By rwe denote the position vector which is further expressed by spherical coordinates r= (r,Ω) = (r,ϑ,ϕ ). We approximate the Earth by a sphere Gwith surface ∂Gof radiusa. We assume that the sphere Ghas constant magnetic permeability µ0and positive, finite, continuous 3-D varying electrical conductivity σ(r)∈L+ ∞(G)∩C1(G) (see Appendix F.2 for a summary of used functional spaces). It is also useful to introduce the electrical resistivity ρ(r) = 1/σ(r), ρ(r)∈L+ ∞(G)∩C1(G). The magnetic field in Gis induced by time-varying electrical currents j(e)in the ionosphere and the magnetosphere which are separated from the Earth by a perfectly insulating spherical atmosphere A. The model configuration is sketched in Figure 2.1. The characteristic time scales of the external current variations range from several hours to hundreds of days, so the quasi-static approximation of Maxwell’s equations is applicable. Under these assumptions, the equation of electromagnetic induction can be derived directly from Maxwell’s equations (Parkinson & Hutton, 1989): µ0∂B ∂t+ curl (ρcurlB) = 0 inG, (2.1) where B(r;t)∈C2(G)3×C1(/angbracketleft0,∞))3is the vector of magnetic induction. Equation (2.1) is a parabolic differential equation. The initial value and the boundary value conditions ofBhave to be specified. The initial value0Bat timet= 0 is specified by a function fromC2(G)3that satisfies the divergence-free condition div0B= 0 inG. (2.2) Applying the divergence operator on equation (2.1) yields ∂(divB) ∂t= 0 inG. (2.3) Equation (2.3) and the choice of0Bin accordance with (2.2) assure that div B= 0 at any timet>0, i.e. the divergence-free condition on Bis implicitly satisfied by the equation of electromagnetic induction (2.1). 10 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH A σ = 0 δG r ρ(r,ϑ,ϕ)Gϑ y xz a ϕr j(e)n Figure 2.1: A schematic figure of model configuration. Heterogeneous conductive sphere Gof radius a with surface ∂Gis surrounded by an insulating spherical layer A. Above A, a system of primary inducing electrical currents j(e)is prescribed. The general form of the boundary conditions imposed on the magnetic field by Maxwell’s equations on any material interface is, [n·B]+ −= 0, (2.4) [n×H]+ −=jS, (2.5) where H,jS, andnare the magnetic intensity vector, the density of free surface currents, and the unit normal vector on the interface, respectively. Because of the assumption of constant magnetic permeability µ0in the material relation B=µ0H, (2.6) 2.1. FORMULATION OF THE INITIAL-VALUE PROBLEM 11 and because no free surface currents are generated on the surfaces of finitely conducting, or insulating bodies, the boundary conditions imposed on Bon∂Greduce to [B]+ −= 0. (2.7) The magnetic induction vector Bchanges continuously across the spherical surface ∂G. The magnetic induction in the insulating atmosphere Ais a divergence-free and irro- tational field, and can be represented at any time t≥0 by a gradient of scalar magnetic potentialU(r;t), B(r;t) =−gradU(r;t) inA, (2.8) whereUsatisfies the Laplace equation, ∇2U(r;t) = 0 inA. (2.9) Expanding the magnetic potential into a series of spherical harmonics Yjm(Ω) (Var- shalovich et al., 1989) allows us to write the analytical solution of (2.9) as a sum of two parts, U=U(e)+U(i). (2.10) The external part, U(e)(r;t) =a∞/summationdisplay j=1j/summationdisplay m=−jG(e) jm(t)/parenleftBigr a/parenrightBigj Yjm(Ω), (2.11) represents the magnetic field induced by ionospheric and magnetospheric currents j(e) above the atmospheric layer A. The internal part, U(i)(r;t) =a∞/summationdisplay j=1j/summationdisplay m=−jG(i) jm(t)/parenleftBiga r/parenrightBigj+1 Yjm(Ω), (2.12) represents the magnetic field of induced currents inside the sphere G(e.g., Langel, 1987). Substituting from (2.10), (2.11), and (2.12) into (2.8), and taking into account the formula for gradient of spherical harmonic functions (A 15), yields B(r;t) =B(e)(r;t) +B(i)(r;t) inA, (2.13) where B(e)(r;t) = −/summationdisplay jmG(e) jm(t)/parenleftBigr a/parenrightBigj−1/bracketleftBig jS(−1) jm(Ω) + S(1) jm(Ω)/bracketrightBig , (2.14) B(i)(r;t) =/summationdisplay jmG(i) jm(t)/parenleftBiga r/parenrightBigj+2/bracketleftBig (j+ 1)S(−1) jm(Ω)−S(1) jm(Ω)/bracketrightBig . (2.15) 12 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH The vector spherical harmonic functions S(λ) jm(Ω),λ=−1,0,1, are introduced in Ap- pendix A. The limits in the sums over jandmremain the same as in equation (2.11) and will be omitted from now on. The downward, northward, and eastward components of magnetic induction in the atmosphere can be derived from the definitions (A 4–A 6), Z(r;t) = −er·B(r;t) = =/summationdisplay jm/bracketleftbigg jG(e) jm(t)/parenleftBigr a/parenrightBigj−1 −(j+ 1)G(i) jm(t)/parenleftBiga r/parenrightBigj+2/bracketrightbigg Yjm(Ω),(2.16) X(r;t) = −eϑ·B(r;t) = =/summationdisplay jm/bracketleftbigg G(e) jm(t)/parenleftBigr a/parenrightBigj−1 +G(i) jm(t)/parenleftBiga r/parenrightBigj+2/bracketrightbigg∂Yjm ∂ϑ(Ω), (2.17) Y(r;t) = eϕ·B(r;t) = =−1 sinϑ/summationdisplay jm/bracketleftbigg G(e) jm(t)/parenleftBigr a/parenrightBigj−1 +G(i) jm(t)/parenleftBiga r/parenrightBigj+2/bracketrightbigg∂Yjm ∂ϕ(Ω). (2.18) The unit vectors corresponding to the radial coordinate r, colatitude ϑ, and longitude ϕ, are denoted by er,eϑ, andeϕ. Spherical harmonic coefficients of the external and internal potential,G(e) jm(t) andG(i) jm(t), can be determined from components X,Y,Zmeasured at ground based observatories or at a low-altitude satellite orbit. This problem, although important in 3-D inversions for the conductivity structure, is not discussed in this work. To express the boundary condition at the Earth’s surface, we expand the magnetic induction vector inside the sphere Ginto a series of vector spherical harmonics, B(r,Ω;t) =∞/summationdisplay j=1j/summationdisplay m=−j1/summationdisplay λ=−1B(λ) jm(r;t)S(λ) jm(Ω), (2.19) whereB(λ) jm(r;t) are complex functions of radius and time. Because we use complex pa- rameterization to represent real functions, the property (A 11) of the spherical harmonics yields B(λ) j−m(r;t) = (−1)mB(λ) jm(r;t). (2.20) Assuming the continuity of Bforr=ain equations (2.13–2.15) and (2.19) leads to B(0) jm(a;t) = 0, (2.21) B(−1) jm(a;t) = −/bracketleftBig jG(e) jm(t)−(j+ 1)G(i) jm(t)/bracketrightBig , (2.22) B(1) jm(a;t) = −/bracketleftBig G(e) jm(t) +G(i) jm(t)/bracketrightBig . (2.23) 2.1. FORMULATION OF THE INITIAL-VALUE PROBLEM 13 The spherical harmonic coefficients appropriate to the toroidal part, B(0) jm, vanish at the boundary∂G, the external and internal parts of the magnetic potential are present in the vertical spheroidal components B(−1) jmand in the horizontal spheroidal components B(1) jm. By eliminating G(i) jm(t) from (2.22) and (2.23), we can reduce the boundary conditions to B(0) jm(a;t) = 0, (2.24) B(−1) jm(a;t) + (j+ 1)B(1) jm(a;t) = −(2j+ 1)G(e) jm(t). (2.25) The boundary conditions (2.24–2.25) are connected to the spherical harmonic expansion ofB, and cannot be expressed as simple Dirichlet or Neumann boundary conditions imposed on B. Therefore, they are difficult to implement in other numerical schemes, such as 3-D finite elements, or 3-D finite differences. In such cases, the Laplace equation (2.9) is solved numerically, and boundary condition on potential Uis prescribed on the outer surface ofA. This approach is implemented in the 3-D finite element approach in Chapter 3, see also Everett & Schultz (1996) and Uyeshima & Schultz (2000). However, using the orthogonality of the spherical harmonic functions (A 8), it is possible to reformulate equations (2.24–2.25) as a set of integral equations /integraldisplay ΩB(a,Ω;t)·S(0) jm(Ω)dΩ = 0, (2.26) /integraldisplay ΩB(a,Ω;t)·/bracketleftbigg S(−1) jm(Ω) +1 jS(1) jm(Ω)/bracketrightbigg dΩ = −(2j+ 1)G(e) jm(t) (2.27) which hold for each j≥1 andm, such that |m| ≤j. Bar denotes complex conjugation. Now we can finally summarize the formulation of the electromagnetic induction prob- lem in the time-domain in the classical sense: Letµ0>0 andρ(r)∈L+ ∞(G)∩C1(G). Let the initial value0B(r)∈C2(G)3be a divergence-free function, div0B= 0 inG. LetU(e)(r,Ω;t) be the external part of magnetic potential in the atmosphere defined by (2.11). Find such a function B(r;t)∈C2(G)3×C1(/angbracketleft0,∞))3that satisfies µ0∂B ∂t+ curl (ρcurlB) = 0 in G, (2.28) B(r; 0) =0B(r) inG, (2.29) (2.30) 14 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH and∀j≥1,|m| ≤j, /integraldisplay ΩB(a,Ω;t)·S(0) jm(Ω)dΩ = 0, (2.31) /integraldisplay ΩB(a,Ω;t)·/bracketleftbigg S(−1) jm(Ω) +1 jS(1) jm(Ω)/bracketrightbigg dΩ = −(2j+ 1)G(e) jm(t). (2.32) It is also possible to extend the classical formulation (2.28–2.32) to the case of dis- continuous resistivity in G. Let us assume internal surfaces Γ i, where resistivity ρ(r) changes discontinuously, separating the sphere Ginto subdomains Gi, whereρ(r) is con- tinuous. Then the equation (2.28) holds in each subdomain Gi, and continuity of B(2.7) is imposed on the internal interfaces Γ i. 2.1.2 Integral (weak) formulation To reformulate the problem of electromagnetic induction in a heterogeneous sphere in the integral sense, we introduce the solution spaces Hcurl,Hcurl,0(Kˇ r´ ıˇ zek & Neittaanm¨ aki, 1990), and the L2scalar product ( ·,·) (see Appendix F.2). Then we can pose the problem as follows: Letµ0>0 andρ(r)∈L+ ∞(G). Let the initial value0B(r)∈Hcurlbe a divergence-free function, div0B= 0 inG. LetU(e)(r,Ω;t) be the external part of magnetic potential in the atmosphere defined by (2.11). Find B(r;t)∈Hcurl×C1(/angbracketleft0,∞))3, such that µ0/parenleftbigg∂B ∂t,δB/parenrightbigg +a(B,δB) = 0 ∀δB∈Hcurl,0, (2.33) B(r; 0) =0B(r) inG, (2.34) and∀j≥1,|m| ≤j, /integraldisplay ΩB(a,Ω;t)·S(0) jm(Ω)dΩ = 0, (2.35) /integraldisplay ΩB(a,Ω;t)·/bracketleftbigg S(−1) jm(Ω) +1 jS(1) jm(Ω)/bracketrightbigg dΩ = −(2j+ 1)G(e) jm(t). (2.36) 2.1. FORMULATION OF THE INITIAL-VALUE PROBLEM 15 The sesquilinear form a(·,·) is defined as a(B,δB) =/integraldisplay Gρ(curlB)·(curlδB)dV. (2.37) Let us show that if there exists a weak solution of (2.33–2.36), such that B∈C2(G)3× C1(/angbracketleft0,∞))3, forρ∈C1(G) and0B∈C2(G)3, it is also the classical solution of (2.28– 2.32). Since C2(G)3⊂Hcurl, we can take δB∈C2(G), such that er×δB= 0 on∂G, and apply on (2.33) Green’s theorem, /integraldisplay Gcurl (hcurlf)·gdV=/integraldisplay Gh(curlf)·(curlg)dV−/integraldisplay ∂Gh(curlf)·(er×g)dS,(2.38) valid for any h∈C1(G),f,g∈C2(G)3. The surface integral is canceled because of the choice ofδB, and we are left with /integraldisplay G/bracketleftbigg µ0∂B ∂t+ curl (ρcurlB)/bracketrightbigg ·δBdV= 0. (2.39) This can be satisfied for all δB∈C2(G),er×δB= 0 on∂G, only if µ0∂B ∂t+ curl (ρcurlB) = 0 inG. (2.40) The integral formulation also implicitly satisfies the divergence-free constraint on B. Let us introduce an auxiliary function ξ∈D0(G), and construct a particular test func- tionsδB= gradξ. Note that δB∈Hcurl,0because grad ξ∈L2(G)3, curl gradξ= 0, and er×gradξ= 0 on∂G. Then equation (2.33) yields /integraldisplay G∂B ∂t·gradξdV = 0 ∀ξ∈D0(G). (2.41) BecauseHcurl⊂L2(G)3andD0(G)⊂L2(G), we can use another Green’s theorem which holds for any square integrable scalar hand vector f, /integraldisplay Gf·gradhdV =−/integraldisplay G(divf)hdV +/integraldisplay ∂G(er·f)hdS. (2.42) 16 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH The surface integral is zero because ξ= 0 on∂G, and the time and space derivatives are interchangeable, therefore we can write /integraldisplay G∂(divB) ∂tξdV = 0 ∀ξ∈D0(G). (2.43) This can be satisfied only with ∂(divB) ∂t= 0 inG. (2.44) However, numerical errors can introduce spurious magnetic monopoles during the time integration. This can be avoided either by adding the divergence-free constraint to (2.33) in the form of the Lagrange multipliers (Martinec, 1999a), or by a projection to the divergence-free solution (Uyeshima & Schultz, 2000). Although a modified formulation of the presented method which uses the former approach to assure solenoidality of Bwas also developed, we have not observed any recognizable deviations from the divergence- free solution while testing the code for a wide selection of resistivity models and external source currents, and the constraint was dropped from the algorithm in favour of increased speed and lower memory requirements of the numerical method. 2.1.3 Time integration scheme Implicit time integration schemes for solving parabolic differential equations are gener- ally preferred to explicit ones (Press et al., 1992). However, with fully 3-D resistivity distribution, and due to the traditional spherical harmonic parameterization used in the angular coordinates, the form a(·,·) leads to a very large matrix. Martinec (1999a) pro- posed a method for a fast iterative solution of such a system, but it would still require to compute many iterations at each level of the time discretization. We present a different, semi-implicit approach, based on splitting of the form a(·,·) into two parts. Let us express the resistivity as a sum of two components, ρ(r,Ω) =ρ0(r) +ρ1(r,Ω), (2.45) whereρ0(r)∈L+ ∞(G) is a spherically symmetric, finite, and positive resistivity model, andρ1(r,Ω) represents the deviations (not necessarily small) of the actual resistivity ρ from the spherically symmetric model. We experienced the best numerical results by prescribing ρ0as a maximum of ρover angular coordinates, ρ0(r) = max Ωρ(r,Ω). (2.46) 2.1. FORMULATION OF THE INITIAL-VALUE PROBLEM 17 Now we can split the sesquilinear form (2.37) into two parts corresponding to ρ0andρ1, respectively, a(B,δB) =a0(B,δB) +a1(B,δB), (2.47) where a0(B,δB) =/integraldisplay Gρ0(curlB)·(curlδB)dV, (2.48) a1(B,δB) =/integraldisplay Gρ1(curlB)·(curlδB)dV. (2.49) As we will show later, form a0(·,·) is represented by a banded matrix with 9 non-zero diagonals. Therefore, it is easy to treat it implicitly in the time integration scheme, while more complex form a1(·,·) is considered on the right hand side of equation (2.33), and is evaluated using the solution Bfrom the previous time step. Approximating the time derivative by the difference between the values of magnetic field at two following time steps, ∂B ∂t≈i+1B−iB ∆t, (2.50) whereiB=B(r;ti)∈Hcurl, and the time step ∆ t=ti+1−tiis considered constant, we can formulate the time integration scheme as a sequence of elliptic problems with integral boundary conditions: Findi+1B∈Hcurl∀i= 0,1,..., such that µ0 ∆t/parenleftbigi+1B,δB/parenrightbig +a0(i+1B,δB) =µ0 ∆t/parenleftbigiB,δB/parenrightbig −a1(iB,δB)∀δB∈Hcurl,0,(2.51) and∀j≥1,|m| ≤j, /integraldisplay Ωi+1B(a,Ω)·S(0) jm(Ω)dΩ = 0, (2.52) /integraldisplay Ωi+1B(a,Ω)·/bracketleftbigg S(−1) jm(Ω) +1 jS(1) jm(Ω)/bracketrightbigg dΩ = −(2j+ 1)G(e) jm(ti+1).(2.53) The initial value0Bfrom (2.34) is used on the right hand side of (2.51) for i= 0. Note that the solutioni+1Band the test functions δBare from different functional spaces Hcurl andHcurl,0, respectively. The existence and uniqueness of the elliptic problem (2.51–2.53) are discussed in Appendix B. 18 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH ψ11 0ψ2 r2 k−1rkrk+1rPPr+1= a r r1= 0ψk−1ψkψk+1ψPψP+1 Figure 2.2: One-dimensional piecewise linear finite element functions. 2.1.4 Spherical harmonic-finite element parameterization The expansion of the magnetic induction Binto a series of vector spherical harmonics has been introduced by equation (2.19). To complete the spatial discretization of the problem (2.51–2.53), we need to define a parameterization over the radial coordinate. We use piecewise linear finite elements for this purpose. Let us divide the interval /angbracketleft0,a/angbracketright intoPsubintervals Ik=/angbracketleftrk,rk+1/angbracketright, where 0 = r1< r 2< ... < r P< rP+1=a. In each interval Ik, there are only two non-zero finite element functions, namely ψk(r) =rk+1−r hk, ψ k+1(r) =r−rk hk, (2.54) wherehk=rk+1−rkis the length of the interval Ik(see Figure 2.2). Now we are ready to introduce the discrete approximation of the functional space Hcurl, Hh curl=/braceleftBigg B/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleB=jmax/summationdisplay j=1j/summationdisplay m=−jP+1/summationdisplay k=11/summationdisplay λ=−1B(λ) jm,kψkS(λ) jm/bracerightBigg . (2.55) The spherical harmonic expansion is truncated at finite degree jmax. Note that the spherical harmonic coefficients B(λ) jm,kmust satisfy the discrete form of (2.20), B(λ) j−m,k= (−1)mB(λ) jm,k, (2.56) since Bis real. By excluding the basis functions ψP+1S(0) jmandψP+1S(1) jmwhich have non-zero horizontal components at the surface we get an approximation of the space 2.1. FORMULATION OF THE INITIAL-VALUE PROBLEM 19 Hcurl,0, Hh curl,0=/braceleftBigg δB/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingleδB=/summationdisplay jm/bracketleftBiggP/summationdisplay k=11/summationdisplay λ=−1δB(λ) jm,kψkS(λ) jm+δB(−1) jm,P +1ψP+1S(−1) jm/bracketrightBigg/bracerightBigg .(2.57) The Galerkin approximation of the solution of (2.51–2.53) is then such a functioni+1B∈ Hh curlthat satisfies (2.51) for all test functions δB∈Hh curl,0, and the boundary conditions (2.52–2.53) for j= 1,...,j max,m= 0,...,j . Now we aim to derive the spherical harmonic-finite element approximations of the L2 scalar product ( ·,·), and of the sesquilinear forms a0(·,·) anda1(·,·). Since Bis a real function, we can replace it by complex conjugation and make use of the orthogonality of the spherical harmonics (A 8). Separating the radial and angular parts of the volume integral,/integraltext GdV=/integraltexta 0r2dr/integraltext ΩdΩ, and substituting for BandδBthe basis functions ψk(r)S(λ) jm(Ω), yields for the L2scalar product, /parenleftBig ψkS(λ) jm,ψk/primeS(λ/prime) j/primem/prime/parenrightBig =δjj/primeδmm/primeδλλ/primeNjλ/bracketleftBig I(2) kδk(k/prime−1)+/parenleftBig I(1) k+I(3) k−1/parenrightBig δkk/prime+I(2) kδk(k/prime+1)/bracketrightBig . (2.58) The normNjλis given by (A 9). In the previous equations we have introduced the inte- grals of the finite element products over intervals Ikwhich can be expressed analytically, I(1) k=rk+1/integraldisplay rkψk(r)ψk(r)r2dr=hk 30/parenleftbig r2 k+1+ 3rkrk+1+ 6r2 k/parenrightbig , (2.59) I(2) k=rk+1/integraldisplay rkψk(r)ψk+1(r)r2dr=hk 60/parenleftbig 3r2 k+1+ 4rkrk+1+ 3r2 k/parenrightbig , (2.60) I(3) k=rk+1/integraldisplay rkψk+1(r)ψk+1(r)r2dr=hk 30/parenleftbig 6r2 k+1+ 3rkrk+1+r2 k/parenrightbig . (2.61) These are defined for k= 1,...,P . To keep unified notation of formula (2.58) even for the lowermost ( k,k/prime= 1) and for the uppermost ( k,k/prime=P) layers, we formally define I(·) 0=I(·) P+1= 0. Before we express the approximation of the form a0(·,·) we define on each interval Ik functions χk(r) =r/parenleftbiggd dr+1 r/parenrightbigg ψk(r) =rk+1−2r hk, (2.62) χk+1(r) =r/parenleftbiggd dr+1 r/parenrightbigg ψk+1(r) =2r−rk hk. (2.63) 20 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH Then, using (A 16–A 18), we can write the rotation of the basis functions ψkS(λ) jmfrom Hh curlas rcurl/bracketleftBig ψk(r)S(0) jm(Ω)/bracketrightBig =−Πjψk(r)S(−1) jm(Ω)−χk(r)S(1) jm(Ω), (2.64) rcurl/bracketleftBig ψk(r)S(−1) jm(Ω)/bracketrightBig =−ψk(r)S(0) jm(Ω), (2.65) rcurl/bracketleftBig ψk(r)S(1) jm(Ω)/bracketrightBig =χk(r)S(0) jm(Ω). (2.66) Introducing the piecewise constant parameterization of the radial resistivity profile ρ0(r), ρ0(r) =ρ0,k= const for r∈ Ik, (2.67) and substituting (2.64–2.66) into the definition (2.48), yields for the form a0(·,·), a0/parenleftBig ψkS(λ) jm,ψk/primeS(0) j/primem/prime/parenrightBig =δjj/primeδmm/primeδλ0Πjρ0,k/prime/bracketleftbigg/parenleftBig ΠjK(2) k+K(5) k/parenrightBig δk(k/prime−1)+ +/parenleftBig ΠjK(1) k+K(4) k+ Π jK(3) k−1+K(6) k−1/parenrightBig δkk/prime+ +/parenleftBig ΠjK(2) k+K(5) k/parenrightBig δk(k/prime+1)/bracketrightbigg , (2.68) a0/parenleftBig ψkS(λ) jm,ψk/primeS(−1) j/primem/prime/parenrightBig =δjj/primeδmm/primeΠjρ0,k/prime/bracketleftbigg K(2) k−1δk(k/prime−1)δλ−1−K(8) k−1δk(k/prime−1)δλ1+ +/parenleftBig K(1) k+K(3) k−1/parenrightBig δkk/primeδλ−1−/parenleftBig K(7) k+K(10) k−1/parenrightBig δkk/primeδλ1+ +K(2) kδkk/primeδλ−1−K(9) kδkk/primeδλ1/bracketrightbigg , (2.69) a0/parenleftBig ψkS(λ) jm,ψk/primeS(1) j/primem/prime/parenrightBig =δjj/primeδmm/primeΠjρ0,k/prime/bracketleftbigg −K(9) k−1δk(k/prime−1)δλ−1+K(5) k−1δk(k/prime−1)δλ1− −/parenleftBig K(7) k+K(10) k−1/parenrightBig δkk/primeδλ−1+/parenleftBig K(4) k+K(6) k−1/parenrightBig δkk/primeδλ1− −K(8) kδk(k/prime+1)δλ−1+K(5) kδk(k/prime+1)δλ1/bracketrightbigg , (2.70) where Π jis defined in (A 10). We have introduced additional integrals over Ik, namely, K(1) k=rk+1/integraldisplay rkψk(r)ψk(r)dr=hk 3, (2.71) K(2) k=rk+1/integraldisplay rkψk(r)ψk+1(r)dr=hk 6, (2.72) 2.1. FORMULATION OF THE INITIAL-VALUE PROBLEM 21 K(3) k=rk+1/integraldisplay rkψk+1(r)ψk+1(r)dr=hk 3, (2.73) K(4) k=rk+1/integraldisplay rkχk(r)χk(r)dr=1 3hk/parenleftbig r2 k+1−2rk+1rk+ 4r2 k/parenrightbig , (2.74) K(5) k=rk+1/integraldisplay rkχk(r)χk+1(r)dr=−1 3hk/parenleftbig r2 k+1+rk+1rk+r2 k/parenrightbig , (2.75) K(6) k=rk+1/integraldisplay rkχk+1(r)χk+1(r)dr=1 3hk/parenleftbig 4r2 k+1−2rk+1rk+r2 k/parenrightbig , (2.76) K(7) k=rk+1/integraldisplay rkψk(r)χk(r)dr=1 6(rk+1−4rk), (2.77) K(8) k=rk+1/integraldisplay rkψk+1(r)χk(r)dr=−1 6(rk+1+ 2rk), (2.78) K(9) k=rk+1/integraldisplay rkψk(r)χk+1(r)dr=1 6(2rk+1+rk), (2.79) K(10) k =rk+1/integraldisplay rkψk+1(r)χk+1(r)dr=1 6(4rk+1−rk), (2.80) fork= 1,...,P . Analogously to I(·) k, we formally define K(·) 0=K(·) P+1= 0. Note that both the L2scalar product and the form a0(·,·) are decoupled in indeces jandm, and the coupling in the index kis limited to neighbouring layers. The former is not valid in the approximation of the form a1(·,·) (2.49), where lateral variations of ρ1introduce coupling in jandm. In the evaluation of a1(·,·) we follow generally the technique developed by Martinec (1999a) with minor modifications due to our choice of different system of vector spherical harmonics. We prescribe the resistivity variations on an angular grid of Nϑ×Nϕnodes, and piecewise constant with respect to the radial coordinate, i.e., ρ1(r,ϑp,ϕq) =ρpq 1,k= const for r∈ Ik, (2.81) where colatitudinal grid nodes {ϑp}Nϑ p=1are the roots of Legendre polynomial of degree 22 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH Nϑ, PNϑ(cosϑp) = 0, (2.82) and longitudinal grid nodes {ϕq}Nϕ−1 q=0are equidistant, ϕq=2πq Nϕ. (2.83) Moreover, we require that Nϑ>3jmax/2,Nϕ>3jmax, andNϕis an integer power of 2 (Martinec, 1989). Now we define the following sums on the grid ( ϑp,ϕq), and for all radial indeces k= 1,...,P + 1, Q(1) pq,k =/summationdisplay jmΠjiB(0) jm,kYjm(ϑp,ϕq), (2.84) Q(2) pq,k =/summationdisplay jmiB(0) jm,k∂Yjm ∂ϑ(ϑp,ϕq), (2.85) Q(3) pq,k =/summationdisplay jmiB(0) jm,k1 sinϑp∂Yjm ∂ϕ(ϑp,ϕq), (2.86) Q(4) pq,k =/summationdisplay jmiB(−1) jm,k∂Yjm ∂ϑ(ϑp,ϕq), (2.87) Q(5) pq,k =/summationdisplay jmiB(−1) jm,k1 sinϑp∂Yjm ∂ϕ(ϑp,ϕq), (2.88) Q(6) pq,k =/summationdisplay jmiB(1) jm,k∂Yjm ∂ϑ(ϑp,ϕq), (2.89) Q(7) pq,k =/summationdisplay jmiB(1) jm,k1 sinϑp∂Yjm ∂ϕ(ϑp,ϕq). (2.90) Because of the exponential function exp ( −imϕ q) in the spherical harmonic functions Yjm, the summation over mfor each node ϕqcan be computed by the FFT algorithm (Martinec, 1989; Press et al., 1992). The sumsQ(α) pq,kare then multiplied by discretized resistivity and spherical harmonics, or their angular derivatives, and numerically integrated over the angular coordinates, C(11) jm,k =2π NϕΠjNϑ/summationdisplay p=1Nϕ−1/summationdisplay q=0wpρpq 1,kQ(1) pq,kYjm(ϑp,ϕq), (2.91) 2.1. FORMULATION OF THE INITIAL-VALUE PROBLEM 23 C(α2) jm,k =2π NϕNϑ/summationdisplay p=1Nϕ−1/summationdisplay q=0wpρpq 1,kQ(α) pq,k∂Yjm ∂ϑ(ϑp,ϕq), (2.92) C(α3) jm,k =2π NϕNϑ/summationdisplay p=1Nϕ−1/summationdisplay q=0wpρpq 1,kQ(α) pq,k1 sinϑp∂Yjm ∂ϕ(ϑp,ϕq), (2.93) D(11) jm,k =2π NϕΠjNϑ/summationdisplay p=1Nϕ−1/summationdisplay q=0wpρpq 1,kQ(1) pq,k+1Yjm(ϑp,ϕq), (2.94) D(α2) jm,k =2π NϕNϑ/summationdisplay p=1Nϕ−1/summationdisplay q=0wpρpq 1,kQ(α) pq,k+1∂Yjm ∂ϑ(ϑp,ϕq), (2.95) D(α3) jm,k =2π NϕNϑ/summationdisplay p=1Nϕ−1/summationdisplay q=0wpρpq 1,kQ(α) pq,k+11 sinϑp∂Yjm ∂ϕ(ϑp,ϕq), (2.96) for allj= 1,...,j max,m= 0,...,j , andk= 1,...,P . Indexαranges from 2 to 7. In the previous equations we employed the Gauss-Legendre formula with weights wp(Press et al., 1992) for numerical integration over ϑ, π/integraldisplay 0f(ϑ) sinϑdϑ≈Nϑ/summationdisplay p=1wpf(ϑp), (2.97) and the discrete FFT over ϕ, 2π/integraldisplay 0f(ϕ)eimϕdϕ≈2π NϕNϕ−1/summationdisplay q=1f(ϕq)eimϕ q. (2.98) The spherical harmonic-finite element approximation of the form a1(·,·) then can be evaluated for each basis function from Hh curl,0by combining the angular integrals (2.91– 2.96) with the radial integrals (2.71–2.80), and taking into account the rotation formulae (2.64–2.66), a1(iB,ψkS(0) jm) =/bracketleftBig C(11) jm,k−1K(2) k−1+D(11) jm,k−1K(3) k−1+C(11) jm,kK(1) k+D(11) jm,kK(2) k/bracketrightBig + +/bracketleftBig C(22) jm,k−1K(5) k−1+D(22) jm,k−1K(6) k−1+C(22) jm,kK(4) k+D(22) jm,kK(5) k/bracketrightBig + +/bracketleftBig C(33) jm,k−1K(5) k−1+D(33) jm,k−1K(6) k−1+C(33) jm,kK(4) k+D(33) jm,kK(5) k/bracketrightBig + +/bracketleftBig C(43) jm,k−1K(9) k−1+D(43) jm,k−1K(10) k−1+C(43) jm,kK(7) k+D(43) jm,kK(8) k/bracketrightBig − 24 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH −/bracketleftBig C(52) jm,k−1K(9) k−1+D(52) jm,k−1K(10) k−1+C(52) jm,kK(7) k+D(52) jm,kK(8) k/bracketrightBig − −/bracketleftBig C(63) jm,k−1K(5) k−1+D(63) jm,k−1K(6) k−1+C(63) jm,kK(4) k+D(63) jm,kK(5) k/bracketrightBig + +/bracketleftBig C(72) jm,k−1K(5) k−1+D(72) jm,k−1K(6) k−1+C(72) jm,kK(4) k+D(72) jm,kK(5) k/bracketrightBig ,(2.99) a1(iB,ψkS(−1) jm) =/bracketleftBig C(32) jm,k−1K(8) k−1+D(32) jm,k−1K(10) k−1+C(32) jm,kK(7) k+D(32) jm,kK(9) k/bracketrightBig − −/bracketleftBig C(23) jm,k−1K(8) k−1+D(23) jm,k−1K(10) k−1+C(23) jm,kK(7) k+D(23) jm,kK(9) k/bracketrightBig + +/bracketleftBig C(42) jm,k−1K(2) k−1+D(42) jm,k−1K(3) k−1+C(42) jm,kK(1) k+D(42) jm,kK(2) k/bracketrightBig + +/bracketleftBig C(53) jm,k−1K(2) k−1+D(53) jm,k−1K(3) k−1+C(53) jm,kK(1) k+D(53) jm,kK(2) k/bracketrightBig − −/bracketleftBig C(62) jm,k−1K(8) k−1+D(62) jm,k−1K(10) k−1+C(62) jm,kK(7) k+D(62) jm,kK(9) k/bracketrightBig − −/bracketleftBig C(73) jm,k−1K(8) k−1+D(73) jm,k−1K(10) k−1+C(73) jm,kK(7) k+D(73) jm,kK(9) k/bracketrightBig ,(2.100) a1(iB,ψkS(1) jm) =/bracketleftBig C(23) jm,k−1K(5) k−1+D(23) jm,k−1K(6) k−1+C(23) jm,kK(4) k+D(23) jm,kK(5) k/bracketrightBig − −/bracketleftBig C(32) jm,k−1K(5) k−1+D(32) jm,k−1K(6) k−1+C(32) jm,kK(4) k+D(32) jm,kK(5) k/bracketrightBig − −/bracketleftBig C(42) jm,k−1K(9) k−1+D(42) jm,k−1K(10) k−1+C(42) jm,kK(7) k+D(42) jm,kK(8) k/bracketrightBig − −/bracketleftBig C(53) jm,k−1K(9) k−1+D(53) jm,k−1K(10) k−1+C(53) jm,kK(7) k+D(53) jm,kK(8) k/bracketrightBig + +/bracketleftBig C(62) jm,k−1K(5) k−1+D(62) jm,k−1K(6) k−1+C(62) jm,kK(4) k+D(62) jm,kK(5) k/bracketrightBig + +/bracketleftBig C(73) jm,k−1K(5) k−1+D(73) jm,k−1K(6) k−1+C(73) jm,kK(4) k+D(73) jm,kK(5) k/bracketrightBig .(2.101) Let us order the unknown coefficientsi+1B(λ) jm,kof the magnetic field at the ( i+ 1)-th time step into a complex column vector xT=/braceleftBigg/bracketleftbigg/parenleftBig i+1B(0) jm,k,i+1B(−1) jm,k,i+1B(1) jm,k/parenrightBigP+1 k=1/bracketrightbiggj m=0/bracerightBiggjmax j=1, (2.102) with dimension d= 3(P+ 1)jmax(jmax+ 3)/2. The complex coefficients for negative mare not included in vector x, but, when needed in the evaluation of a1(·,·), they are computed using the symmetry (2.56). Now we construct the Galerkin system of dlinear complex equations which represents a discrete approximation of (2.51–2.53), A ·x=y. (2.103) 2.2. VALIDATION 25 The elements of matrix Aand the right-hand side vector yare defined as an/primen=µ0 ∆t/parenleftBig ψkS(λ) jm,ψk/primeS(λ/prime) j/primem/prime/parenrightBig +a0/parenleftBig ψkS(λ) jm,ψk/primeS(λ/prime) j/primem/prime/parenrightBig , (2.104) yn/prime=µ0 ∆tjmax/summationdisplay j=1j/summationdisplay m=−jP+1/summationdisplay k=11/summationdisplay λ=−1iB(λ) jm,k/parenleftBig ψkS(λ) jm,ψk/primeS(λ/prime) j/primem/prime/parenrightBig −a1/parenleftBig iB,ψk/primeS(λ/prime) j/primem/prime/parenrightBig ,(2.105) where the multi-indeces nandn/primeare considered from intervals n={j,m,k,λ }={{1,...,j max},{0,...,j },{1,...,P + 1},{0,−1,1}}, (2.106) n/prime={j/prime,m/prime,k/prime,λ/prime}={{1,...,j max},{0,...,j/prime},{1,...,P +δλ/prime−1},{0,−1,1}}.(2.107) Because the test functions ψP+1S(0) jmandψP+1S(1) jmwere excluded from Hh curl,0, there are two missing equations for each j/primem/primein equations (2.104–2.107). Hence, the system must be completed by the boundary conditions (2.24–2.25). The constraint on the toroidal field on the surface ∂Gis imposed by prescribing an/primen=δjj/primeδmm/primeδk(P+1)δλ0, y n/prime= 0, (2.108) forn/prime={{1,...,j max},{0,...,j/prime},P+ 1,0}. The load by external magnetic potential is implemented by setting an/primen=δjj/primeδmm/primeδk(P+1)[(δλ−1+ (j+ 1)δλ1], y n/prime=−(2j/prime+ 1)G(e) j/primem/prime(ti+1), (2.109) forn/prime={{1,...,j max},{0,...,j/prime},P+ 1,1}. Matrix Ais a complex 9-banded matrix which depends only on the length of the time step ∆ tand on the radial resistivity structure ρ0. Since the linear system (2.103) is repeatedly solved in each time step, it is effective to use the LU-decomposition of A. Using a subroutine provided by LAPACK (Anderson et al., 1999), for jmax= 40 and P= 100, it takes only about 10 s on a 500 MHz PC to decompose A. Evaluation of the right-hand side yis based on (2.84–2.96) and (2.99–2.101), and uses the solutioniBfrom the previous time step. It is the most time-consuming part of the algorithm, and takes about 20 s for each time step at given spatial resolution. 2.2 Validation 2.2.1 Nested-sphere conductivity models In this section we validate the presented time-domain spectral-finite element approach using the conductivity models consisting of a homogeneous sphere with a nested spherical inclusion of different conductivity. The configuration of the nested-sphere conductivity 26 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH σ1σ2ϑd b da Figure 2.3: Nested-sphere conductivity model (cross-section). A homogeneous sphere with radius b and conductivity σ2is nested in the host sphere with radius aand conductivity σ1. The centre of the inclusion is positioned at a distance dfrom the centre of the host sphere at colatitude ϑdand longitude ϕd. model is shown in Figure 2.3. This class of models with a frequency-domain, semi- analytical solution developed by Everett & Schultz (1995) for axially symmetric models, and generalized by Martinec (1998) to treat an off-axis spherical inclusion, has become a standard tool for validation of 3-D EM forward solvers (Everett & Schultz, 1996; Martinec, 1999a; Uyeshima & Schultz, 2000). In each homogeneous sphere the equation of EM induction reduces to a complex vector Helmholtz equation, ( ∇2+k2 n)BSA= 0, where k2 n=−iωµ 0σn, and the subscript SAstands for semi-analytical. The magnetic induction vector BSAis expanded into truncated series of eigenfunctions of the Helmholtz operator. The coefficients of the expansion in both spheres are determined from the boundary conditions on the interface by solving a system of linear equations. Although theoretically this method should provide results of almost analytical accu- racy, caution must be paid when using a high-degree spherical harmonic expansion for high frequencies ω, where the method often breaks down numerically. However, this numerical instability has a physical background. As the frequency ωincreases, the pene- tration depth of the damped EM wave decreases, and the solution becomes less sensitive to the position, size, and conductivity of the inclusion. For ω→ ∞ the response of the nested spheres model converges to the response of a homogeneous sphere. In the case of dipolar excitation the response of a homogeneous sphere can be expressed by dipolar terms only. Therefore, the truncation degree jmaxfor a particular frequency band is cho- 2.2. VALIDATION 27 (ω)G10(e)^(t)G10(e) 0100200300400 (nT) 0100200300400 t (h)0100200300400 (nT) 0100200300400 t (h)-2002040(mT s) 1e-091e-081e-081e-071e-061e-051e-04 ω (rad/s)-2002040 1e-091e-081e-081e-071e-061e-051e-04-2002040 1e-091e-081e-081e-071e-061e-051e-04-2002040 1e-091e-081e-081e-071e-061e-051e-04Real Imag Figure 2.4: Analytical model of the geomagnetic storm defined by (2.111) is shown in the left panel. The G(e) 10coefficient derived from the Dst index during a real storm in February 1980 is plotted only for comparison. The Fourier spectrum of the analytical model is shown on the right panel with sampling that is used in the semi-analytical method. Solid and open symbols correspond to truncation degree of 18 and 12, respectively. sen according to two criteria. Firstly, jmaxis small enough to yield a stable solution, i.e., the power of magnetic field is a decreasing function of spherical harmonic degree j. Sec- ondly,jmaxis large enough to parameterize the model sufficiently, i.e., introducing small variations of jmaxin both directions has no significant impact on the obtained solution. We use two following models in the validation. In Model A, the host sphere with radius a= 6371 km has conductivity σ1= 1S/m. The centre of the spherical inclusion with radius b= 3500 km and conductivity σ2= 10S/mis positioned at a distance d= 2700 km from the centre of the host sphere at colatitude ϑd= 40◦and longitude ϕd= 35◦. The Model B is characterized by the same parameters, except for the conductivity of the host sphere which is set to σ1= 0.1S/m, thus increasing the conductivity contrast in the model by two orders of magnitude. 2.2.2 Excitation The presented method allows the excitation by an arbitrary spatio-temporal configuration of the external inducing field. In the following examples we excite the conductivity models by two different signals. In order to compare the results in the frequency-domain we use a harmonic signal 28 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH (Load 1) with simple dipolar spatial structure, G(e) jm(t) =/braceleftBigg A/radicalBig 4π 3sin(ωt) forj= 1, m= 0, 0 otherwise ,(2.110) whereA= 100 nT is the amplitude,/radicalbig 4π/3 is the inverse norm of P10(cosϑ), andω= 3×10−7rad/sis the angular frequency (period T= 242 d). It corresponds to a ring current with periodically varying intensity. The second model (Load 2) of the excitation source used for validation of the time- domain approach is an analytical approximation of the complicated spatio-temporal struc- ture of a geomagnetic storm (Daglis & Kozyra, 2002). We define a time-varying dipolar coefficient of the external magnetic potential as G(e) jm(t) =/braceleftBigg A/radicalBig 4π 3te−α tforj= 1, m= 0, 0 otherwise ,(2.111) whereA= 0.003nT/sis the amplitude, and 1 /α= 48 h is relaxation time of the storm. Regardless of its simplicity, this mathematical model characterizes two basic features of magnetic storms, the P10spatial configuration, and the time evolution with relatively fast onset and slow relaxation (see Figure 2.4). The Fourier spectrum of signal (2.111) has analytical form, ˆG(e) 10(ω) =/radicalbigg 4π 3A 2πα2−ω2−2iαω (α2+ω2)2. (2.112) 2.2.3 Results In the first example we validate the time-domain, spectral-finite element approach against the semi-analytical solution BSAin the frequency-domain. The nested-sphere Model A is excited by the harmonic Load 1 with period T= 242 d. The same configuration was used also by Everett & Schultz (1996) and Martinec (1999a). The long period assures that the EM energy penetrates down to the centre of the sphere, allowing to test the 3-D methods throughout the model. The semi-analytical solution BSAis truncated at jmax= 18. The time-domain solution is started from zero initial condition and evolves for 4 periods, i.e. 968 days. The results are transformed to the frequency-domain using Fourier integration from Tto 4T. In order to avoid the bias from the transient switch-on effect, we exclude the data for t<T . Figure 2.5 shows the convergence of the time integrated solution to the semi-analytical one as the spatial and temporal resolution increases. The response of the nested-sphere model to a slowly varying, long periodic signal is modelled with sufficient accuracy even by the runs with low spatial and temporal resolution. The Bϕcomponent shows the 2.2. VALIDATION 29 highest sensitivity to the spatial parameterization. The time-integration is stable even for 24 h time step, however, the differences from the semi-analytical solution in all three components of Bincrease with ∆ t. In Figure 2.6 we display the time-evolution of the errors with respect to the reference solution BSAas observed on the surface. For the high-resolution model, the initial 100% relative error due to the zero initial state is reduced below 2% after t=T, and goes well below 1% about t= 2Tafter the onset of the time integration. Then the relative error slowly increases, as the numerical errors accumulate. The oscillation of the error is caused by the periodical change of BSA, the local maxima of the error curve correspond to the local minima of the reference harmonic solution. Figure 2.7 shows that applying an additional divergence-free constraint is not neces- sary in the presented formulation. The deviation of the solution from the solenoidal state introduced by numerical errors increases for runs with high lateral resolution ( jmax= 40), decreases for runs with high radial resolution ( P= 100), and is insensitive to the length of time step. Nevertheless, it is negligible for all presented runs. In the case of a high-contrast Model B excited by a simulated geomagnetic storm (Load 2), the reference semi-analytical solution is found for discretely sampled frequencies in the range 10−9–10−4rad/s(see Figure 2.4). The spherical harmonic expansion is truncated atjmax= 18 for frequencies below 5 ×10−6rad/s, and atjmax= 12 for higher frequencies. The time-evolution of BSAis computed by numerical evaluation of Fourier integrals of the frequency-domain solution. Figures 2.8 and 2.9 compare the reference solution with the results of the time-domain, spectral-finite element method. Results obtained by integration with long time steps show considerable delay behind the reference solution. This is a consequence of the explicit treating of the effect of lateral conductivity variations. Although the nested spheres Model B has the same geometry as previous Model A, increased conductivity contrast transfers more energy into the higher degree terms of the spherical harmonic series, and good accuracy is obtained only for the highest truncation degree. 30 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH Br BϑBϕ -40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000Imag RealImagRealImag Real -40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-40-20020(nT) 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-1001020 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000-10010 0200040006000Imag RealImagRealImag Real -40-20020(nT) 0200040006000 r (km)-40-20020(nT) 0200040006000 r (km)-40-20020(nT) 0200040006000 r (km)-40-20020(nT) 0200040006000 r (km)-40-20020(nT) 0200040006000 r (km)-40-20020(nT) 0200040006000 r (km)-40-20020(nT) 0200040006000 r (km)-40-20020(nT) 0200040006000 r (km)-1001020 0200040006000 r (km)-1001020 0200040006000 r (km)-1001020 0200040006000 r (km)-1001020 0200040006000 r (km)-1001020 0200040006000 r (km)-1001020 0200040006000 r (km)-1001020 0200040006000 r (km)-1001020 0200040006000 r (km)-10010 0200040006000 r (km)-10010 0200040006000 r (km)-10010 0200040006000 r (km)-10010 0200040006000 r (km)-10010 0200040006000 r (km)-10010 0200040006000 r (km)-10010 0200040006000 r (km)-10010 0200040006000 r (km)Imag RealImagRealImag Real Figure 2.5: Resolution-dependent convergence of the time-domain method towards the semi-analytical solution for the conductivity Model A excited by Load 1. The components of Bare plotted along radius at fixed colatitude ϑ= 13◦and longitude ϕ= 0◦. The semi-analytical solution is shown by solid line in all plates. Top row: jmax= 40, P= 100, dotted, dash-dotted, and dashed lines correspond to time steps ∆t= 24, 12, and 1 h, respectively. Middle row: P= 100, ∆ t= 1 h, dotted, dash-dotted, and dashed lines correspond to truncation degrees jmax= 10, 20, and 40, respectively. Bottom row: jmax= 40, ∆ t= 1 h, dotted, dash-dotted, and dashed lines correspond to P= 20, 50, and 100 finite elements, respectively. 2.2. VALIDATION 31 100 ||B-BSA||δG/||BSA||δG 0.11110100(%) 0 10000 20000 t (h)100 ||B-BSA||δG/||BSA||δG 0.11110100(%) 0 10000 20000 t (h)100 ||B-BSA||δG/||BSA||δG 0.11110100(%) 0 10000 20000 t (h)100 ||B-BSA||δG/||BSA||δG 0.11110100(%) 0 10000 20000 t (h) Figure 2.6: Relative error of the time-domain solution with respect to the semi-analytical one. The ratio of norms ||B(t)−BSA(t)||L2(∂G)/||BSA(t)||L2(∂G)uses only the values of Bobservable at the surface ∂G. The solid, dashed, dash-dotted, and dotted lines correspond to solutions shown in Figure 2.5 with discretization parameters ( jmax,P, ∆t) set to (40, 100, 1 h), (40, 20, 1 h), (40, 100, 24 h), and (10, 100, 1 h), respectively. ||div B||/||B|| 1e-051e-051e-04(1/m) 0 10000 20000 t (h)||div B||/||B|| 1e-051e-051e-04(1/m) 0 10000 20000 t (h)||div B||/||B|| 1e-051e-051e-04(1/m) 0 10000 20000 t (h)||div B||/||B|| 1e-051e-051e-04(1/m) 0 10000 20000 t (h) Figure 2.7: Deviation from the divergence-free condition introduced by numerical errors during the time integration measured by the ratio of norms ||divB(t)||L2/||B(t)||L2. Same line styles as in Figure 2.6 are used. 32 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH Br BϑBϕ -80-60-40-200(nT) 400050006000-80-60-40-200(nT) 400050006000-80-60-40-200(nT) 400050006000-80-60-40-200(nT) 400050006000-100102030 400050006000-100102030 400050006000-100102030 400050006000-100102030 400050006000-40-30-20-10010 400050006000-40-30-20-10010 400050006000-40-30-20-10010 400050006000-40-30-20-10010 400050006000 -80-60-40-200(nT) 400050006000-80-60-40-200(nT) 400050006000-80-60-40-200(nT) 400050006000-80-60-40-200(nT) 400050006000-100102030 400050006000-100102030 400050006000-100102030 400050006000-100102030 400050006000-40-30-20-10010 400050006000-40-30-20-10010 400050006000-40-30-20-10010 400050006000-40-30-20-10010 400050006000 -80-60-40-200(nT) 400050006000 r (km)-80-60-40-200(nT) 400050006000 r (km)-80-60-40-200(nT) 400050006000 r (km)-80-60-40-200(nT) 400050006000 r (km)-100102030 400050006000 r (km)-100102030 400050006000 r (km)-100102030 400050006000 r (km)-100102030 400050006000 r (km)-40-30-20-10010 400050006000 r (km)-40-30-20-10010 400050006000 r (km)-40-30-20-10010 400050006000 r (km)-40-30-20-10010 400050006000 r (km) Figure 2.8: Resolution-dependent convergence of the time-domain method towards the semi-analytical solution for the conductivity Model B excited by Load 2. The components of Bare plotted at t= 96 h. Colatitude and longitude are fixed at ϑ= 13◦andϕ= 0◦. The semi-analytical solution is shown by solid line in all plates. Top row: jmax= 40, P= 100, dotted, dash-dotted, and dashed lines correspond to time steps ∆ t= 4, 1, and 0 .25 h, respectively. Middle row: P= 100, ∆ t= 0.25 h, dotted, dash-dotted, and dashed lines correspond to truncation degrees jmax= 10, 20, and 40, respectively. Bottom row: jmax= 40, ∆ t= 0.25 h, dotted, dash-dotted, and dashed lines correspond to P= 20, 50, and 100 finite elements, respectively. 2.2. VALIDATION 33 r|B| |B|ϑ|B|ϕ 11010(nT) 020040011010(nT) 020040011010(nT) 020040011010(nT) 02004000.11110 02004000.11110 02004000.11110 02004000.11110 02004000.10.1110 02004000.10.1110 02004000.10.1110 02004000.10.1110 0200400 11010(nT) 020040011010(nT) 020040011010(nT) 020040011010(nT) 02004000.11110 02004000.11110 02004000.11110 02004000.11110 02004000.10.1110 02004000.10.1110 02004000.10.1110 02004000.10.1110 0200400 11010(nT) 0200400 t (h)11010(nT) 0200400 t (h)11010(nT) 0200400 t (h)11010(nT) 0200400 t (h)0.11110 0200400 t (h)0.11110 0200400 t (h)0.11110 0200400 t (h)0.11110 0200400 t (h)0.10.1110 0200400 t (h)0.10.1110 0200400 t (h)0.10.1110 0200400 t (h)0.10.1110 0200400 t (h) Figure 2.9: The same results as in Figure 2.8. The time evolution of the magnetic induction vector is plotted at fixed colatitude ϑ= 13◦and longitude ϕ= 0◦, on the surface of the sphere. 34 2. TIME-DOMAIN, SPECTRAL-FINITE ELEMENT APPROACH 35 Chapter 3 Time-domain, 3-D finite element approach 3.1 A−Φ,Uformulation 3.1.1 Classical formulation We have shown in Chapter 2 that the Galerkin formulation of the EM induction problem in terms of magnetic induction vector Bis based on a discrete approximation of the functional space Hcurl. Although mixed 3-D finite elements conforming to Hcurlcan be constructed (N´ edelec, 1980), it is often advantageous to reformulate the EM induction problems using the scalar and vector electromagnetic potentials which can be parameter- ized by simple nodal finite element functions (Biro & Preis, 1989). Such a technique has already been successfully applied to the electromagnetic induction problem in the Earth’s mantle in the frequency domain (Everett & Schultz, 1996). Here we aim to reformulate it in the time-domain both in the classical and integral sense, and to validate an associated numerical code. The solution domain used in the 3-D finite element method differs from that used in the spectral-finite element approach. Firstly, we reduce the computational requirements by assuming an infinitely conductive core, removing it from the solution domain and replacing it by boundary conditions applied on the core-mantle boundary (CMB). Hence, we approximate the Earth’s mantle by a spherical shell Gwith the inner boundary ∂G1 atr=aCMB and the outer boundary ∂G2atr=a, filled with conductive material of spatially varying conductivity σ(r). Note that it is possible to reformulate the spectral-finite element method from Chap- ter 2 with an infinitely conductive core. However, since the usual choice of homogeneous, highly conductive core introduces no coupling in the spherical harmonic parameteriza- tion, the resulting reduction of the computational burden would be negligible. On the other hand, the 3-D finite element method introduced in this Chapter could incorporate a finitely conductive core, but at the cost of significant increase of computational and memory requirements, since the method does not take advantages of spherically symmet- ric conductivity. Both the finitely and infinitely conductive core approaches are suitable for the realistic configuration of the EM induction problem in the Earth’s mantle. Secondly, since the Laplace equation in the surrounding insulating spherical shell A is solved numerically, not analytically, as in the spectral-finite element method, we define explicitly both surfaces of the outer spherical shell A: the inner boundary ∂A1≡∂G2of 36 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH Gδ2Aδ2 Gδ1Aδ1 y xz r,ϑ,ϕ)G aaO σ( aCMB nnnnσ = 0 A Figure 3.1: A schematic figure of model configuration. Heterogeneous conductive spherical shell G with inner surface ∂G1of radius aCMB and outer surface ∂G2of radius ais surrounded by an insulating spherical layer Awith inner surface ∂A1≡∂G2and outer surface ∂A2of radius aO. Note the orientation of outer normals on the surfaces. radiusr=aand the outer boundary ∂A2of radiusr=aO. Note that the outer unit normal vector nis pointing to and from the Earth’s centre on ∂G1∪∂A1and∂G2∪∂A2, respectively. As can be seen from the analytical formulae (2.14–2.15), the primary inducing field grows asrj−1inA, while the secondary induced field decays as r−j−2, wherejis the spherical harmonic degree. The outer surface ∂A2is placed at such a distance aOthat the field induced by the eddy currents in the conductive mantle is negligible, and only the inducing field due to external current systems can be considered in the boundary condition. A choice of aObetween 4aand 10a(Everett & Schultz, 1996; Uyeshima & Schultz, 2000) yields a decrease in the internal to external field ratio by 2–3 orders of 3.1.A−Φ,UFORMULATION 37 magnitude through the layer A. The vectors of electric intensity and magnetic induction, EandB, in conductors and insulators can be represented by various combinations of scalar and vector potentials. Following Biro & Preis (1989), we use the approach known as A−Φ,Umethod. In this approach, the vector magnetic potential Aand the scalar electrical potential φare used in a conductive material in the forms B(r;t) = −curlA(r;t) inG, (3.1) E(r;t) = −∂A ∂t(r;t)−gradφ(r;t) inG. (3.2) In the time-domain approach it is useful to replace the scalar electrical potential φby its time integral Φ, φ(r;t) =∂ ∂tΦ(r;t) inG, (3.3) since it leads to a symmetric Galerkin matrix. The solenoidal and irrotational magnetic induction Bin the insulator is expressed as a gradient of scalar magnetic potential U, B(r;t) = −gradU(r;t) inA. (3.4) Substitution of (3.1–3.3) and of the Ohm law, j=σE, (3.5) into the quasi static form of the Amp` ere law, curl/parenleftbigg1 µ0B/parenrightbigg =jinG, (3.6) yields curl/parenleftbigg1 µ0curlA/parenrightbigg +σ∂ ∂t(A+ grad Φ) = 0 in G. (3.7) Note that equation (3.7) implicitly satisfies the quasi static continuity equation for the current density, divj=−div/bracketleftbigg σ∂ ∂t(A+ grad Φ)/bracketrightbigg = 0 inG. (3.8) Equation (3.7) does not have a unique solution in terms of Aand Φ. The uniqueness is enforced by applying the Coulomb gauge (Biro & Preis, 1989), divA= 0 inG. (3.9) 38 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH Following Biro & Preis (1989), we apply the Coulomb gauge by addition of a penalty term in equation (3.7), curl/parenleftbigg1 µ0curlA/parenrightbigg −grad/parenleftbigg1 µ0divA/parenrightbigg +σ∂ ∂t(A+ grad Φ) = 0 in G, (3.10) by constraining the solenoidality of current density (3.8) which is not automatically satis- fied by (3.10), and by a homogeneous Dirichlet boundary condition applied on the normal component of the vector magnetic potential on the interface ∂G2, n·A= 0 on∂G2. (3.11) Assuming constant magnetic permeability µ0, equations (3.10) and (3.8) can be further rewritten as, −∇2A+µ0σ∂ ∂t(A+ grad Φ) = 0 in G, (3.12) div/bracketleftbigg µ0σ∂ ∂t(A+ grad Φ)/bracketrightbigg = 0 inG. (3.13) Summed up, we replace the equations (3.7) and (3.9) by equations (3.12–3.13) holding in Gand by the boundary condition (3.11) imposed on ∂G2. Since B= 0 and E= 0 in the infinitely conductive core, we impose the homogeneous Dirichlet boundary conditions on the core-mantle boundary, A= 0 on∂G1, (3.14) Φ = 0 on ∂G1. (3.15) As we have shown in Chapter 2, Bchanges continuously across the Earth’s surface ∂G2≡∂A1. Substitution of (3.1) and (3.4) into (2.7) implies the boundary condition for the vector and scalar magnetic potentials, curlA=−gradUon∂G2≡∂A1. (3.16) Since no electrical current can flow into the insulator, the normal component of the current density must vanish at the Earth’s surface, n·j=n·σE=−σ∂ ∂t(n·A+n·grad Φ) = 0 on ∂G2. (3.17) Taking into account (3.11), the equation (3.17) yields the boundary condition applied on the scalar electric potential only, n·grad Φ = 0 on ∂G2. (3.18) 3.1.A−Φ,UFORMULATION 39 In the insulating region, we solve the Laplace equation for scalar magnetic potential, ∇2U= 0 inA. (3.19) The scalar potential of the primary inducing field of external origin is prescribed at the outer boundary ∂A2, U=U(e)on∂A2. (3.20) To complete the formulation we have to specify the initial conditions. The time evolution of the problem can be started from any initial state (0A,0Φ,0U) that satisfies (3.9) and the boundary conditions (3.11), (3.14–3.16), (3.18), and (3.20). Let us summarize the classical formulation of the A−Φ,Umethod in the time-domain: Letµ0>0 andσ(r)∈L+ ∞(G)∩C1(G). LetU(e)∈C2(∂A2)×C1(/angbracketleft0,∞)) be the magnetic potential of the primary inducing field on ∂A2. Let0A(r)∈C2(G)3, 0Φ(r)∈C2(G), and0U(r)∈C2(A) be the initial values, such that div0A= 0,0A= 0 and0Φ = 0 on∂G1,n·0A=n·grad0Φ = 0 on∂G2, curl0A=−grad0Uon∂G2≡ ∂A1, and0U(aO) =U(e)(Ω; 0). Find such functions A(r;t)∈C2(G)3×C1(/angbracketleft0,∞))3, Φ(r;t)∈C2(G)×C1(/angbracketleft0,∞)), andU(r;t)∈C2(A)×C1(/angbracketleft0,∞)) that satisfy −∇2A+µ0σ∂ ∂t(A+ grad Φ) = 0 in G, (3.21) div/bracketleftbigg µ0σ∂ ∂t(A+ grad Φ)/bracketrightbigg = 0 inG, (3.22) ∇2U= 0 inA, (3.23) A(r; 0) =0AinG, (3.24) Φ(r; 0) =0Φ inG, (3.25) U(r; 0) =0UinA, (3.26) A= 0 on∂G1, (3.27) Φ = 0 on ∂G1, (3.28) curlA=−gradUon∂G2≡∂A1, (3.29) n·A= 0 on∂G2, (3.30) n·grad Φ = 0 on ∂G2, (3.31) U=U(e)on∂A2. (3.32) We refer to Appendix F.2 for definitions of functional spaces. 40 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH 3.1.2 Integral (weak) formulation In order to derive the integral A−Φ,Uformulation of the EM induction problem (3.21– 3.32), we respectively multiply equations (3.21), (3.22), and (3.23) by test functions δA∈W1,2 01(G)3∩C2(G)3,δΦ∈W1,2 01(G)∩C2(G), andδU∈W1,2 02(A)∩C2(A), integrate over domains GandA, and apply the Green theorems /integraldisplay Gf· ∇2gdV =−/integraldisplay Ggradf: (grad g)TdV+/integraldisplay ∂G1∪∂G2f·gradg·ndS− −/integraldisplay ∂G1∪∂G2f·n×curlgdS (3.33) /integraldisplay G(divf)hdV =−/integraldisplay Gf·gradhdV +/integraldisplay ∂G1∪∂G2(n·f)hdS, (3.34) /integraldisplay Af∇2gdV =−/integraldisplay Agradf·gradgdV +/integraldisplay ∂A1∪∂A2fn·gradgdS. (3.35) Integrals over the boundaries ∂G1and∂A2are zero because the homogeneous Dirichlet boundary conditions are imposed on the test functions there, and we obtain /integraldisplay G/bracketleftbigg gradδA: (grad A)T+µ0σδA·∂ ∂t(A+ grad Φ)/bracketrightbigg dV− −/integraldisplay ∂G2δA·gradA·ndS+/integraldisplay ∂G2δA·n×curlAdS= 0,(3.36) /integraldisplay Gµ0σgradδΦ·∂ ∂t(A+ grad Φ)dV−/integraldisplay ∂G2µ0σδΦn·∂ ∂t(A+ grad Φ)dS= 0,(3.37) −/integraldisplay AgradδU·gradU dV +/integraldisplay ∂A1δUn·gradUdS = 0.(3.38) The first surface integral in equation (3.36) and the surface integral in equation (3.37) are zero due to (3.30) and (3.31), respectively. The boundary condition (3.29) is realized by interchanging curl Aand−gradUin the remaining surface integrals in (3.36) and (3.38). These surface integrals couple the solution in the conductive mantle Gto the solution in the insulating layer A(Biro & Preis, 1989, 1990; Everett & Schultz, 1996). By releasing the smoothness constraint on the potentials and test functions, we come to the weak form of the A−Φ,Uformulation: 3.1.A−Φ,UFORMULATION 41 Letµ0>0 andσ(r)∈L+ ∞(G). Let ˆU(r;t) be such a function from W1,2(A)× C1(/angbracketleft0,∞)) that ˆU(aO,Ω;t) =U(e)(Ω;t) on∂A2. Let0A(r)∈W1,2 01(G)3,0Φ(r)∈ W1,2 01(G), and0U(r)∈W1,2(A) be the initial values. Find A(r;t)∈W1,2 01(G)3× C1(/angbracketleft0,∞))3, Φ(r;t)∈W1,2 01(G)×C1(/angbracketleft0,∞)), andU(r;t)∈W1,2(A)×C1(/angbracketleft0,∞)), such thatU−ˆU∈W1,2 02(A) and /integraldisplay G/bracketleftbigg gradδA: (grad A)T+µ0σδA·∂ ∂t(A+ grad Φ)/bracketrightbigg dV− −/integraldisplay ∂G2δA·n×gradU dS = 0∀δA∈W1,2 01(G)3,(3.39) /integraldisplay Gµ0σgradδΦ·∂ ∂t(A+ grad Φ)dV = 0∀δΦ∈W1,2 01(G),(3.40) −/integraldisplay AgradδU·gradU dV −/integraldisplay ∂A1δUn·curlAdS= 0∀δU∈W1,2 02(A).(3.41) 3.1.3 Time integration scheme As we will show in the next section, the use of the 3-D nodal finite elements in the Galerkin discretization of the integral problem (3.39–3.41) leads to a system of linear equations which is much sparser than in the case of the spectral-finite element method introduced in Chapter 2. The reason is that the finite element base functions have local support and introduce less coupling in the Galerkin matrix, than the spherical harmonic functions. Therefore, the application of an implicit time integration method becomes feasible. We use the Crank-Nicolson scheme (Press et al., 1992) which is unconditionally stable and second order accurate in time. Let us denoteiA,iΦ,iUthe solution of (3.39–3.41) at t=ti. The substitution of time derivatives, ∂A ∂t≈i+1A−iA ∆t, (3.42) ∂Φ ∂t≈i+1Φ−iΦ ∆t, (3.43) and of midpoint values, A≈i+1A+iA 2, (3.44) 42 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH Φ≈i+1Φ +iΦ 2, (3.45) U≈i+1U+iU 2, (3.46) into (3.39–3.41) yields the integral formulation discretized in time: Leti+1ˆU(r) be such a function from W1,2(A) thati+1ˆU(aO,Ω) =U(e)(Ω;ti+1) on∂A2. Findi+1A(r)∈W1,2 01(G)3,i+1Φ(r)∈W1,2 01(G), andi+1U(r)∈W1,2(A)∀i= 1,2,..., such thati+1U−i+1ˆU∈W1,2 02(A), and /integraldisplay G/bracketleftbigg1 2gradδA: (gradi+1A)T+µ0σ ∆tδA·/parenleftbigi+1A+ gradi+1Φ/parenrightbig/bracketrightbigg dV−(3.47) −/integraldisplay ∂G21 2δA·n×gradi+1UdS = =/integraldisplay G/bracketleftbigg −1 2gradδA: (gradiA)T+µ0σ ∆tδA·/parenleftbigiA+ gradiΦ/parenrightbig/bracketrightbigg dV + +/integraldisplay ∂G21 2δA·n×gradiUdS ∀δA∈W1,2 01(G)3,(3.48) /integraldisplay Gµ0σ ∆tgradδΦ·/parenleftbigi+1A+ gradi+1Φ/parenrightbig dV = =/integraldisplay Gµ0σ ∆tgradδΦ·/parenleftbigiA+ gradiΦ/parenrightbig dV∀δΦ∈W1,2 01(G),(3.49) −/integraldisplay A1 2gradδU·gradi+1UdV−/integraldisplay ∂A11 2δUn·curli+1AdS= =/integraldisplay A1 2gradδU·gradiUdV +/integraldisplay ∂A11 2δUn·curliAdS∀δU∈W1,2 02(A).(3.50) The existence and uniqueness of the 3-D finite element solution can be proved analo- gously to the spectral-finite element method, as discussed in Appendix B (see also Kˇ r´ ıˇ zek & Neittaanm¨ aki, 1990; Martinec, 1997). Note that the solution does not depend on a particular choice ofi+1ˆU. 3.1.A−Φ,UFORMULATION 43 3.1.4 Discretization using 3-D nodal finite elements Now we will introduce the discrete form of (3.48–3.50) using 3-D piecewise linear nodal finite elements defined on a tetrahedral mesh. A reader not familiar with this technique is advised to read Appendix C that introduces a particular method of tetrahedral mesh generation in a spherical shell, and Appendix D where the 3-D piecewise linear nodal finite elements are defined, before proceeding. Let us assume that the spherical shells GandAare approximated by concave poly- hedra composed respectively of Ni tandNe tnon-overlapping, conforming tetrahedra. The mesh nodes are arranged in such a way that the first Nc nnodes are on the innermost boundary∂G1, nextNi n−Nc n−Ns nnodes are interior nodes of G, followed by Ns nnodes on the interface ∂G2=∂A1,Ne n−Ns n−No ninterior nodes of A, andNo nnodes on the outermost boundary ∂A2. The total number of nodes is then Ni n+Ne n−Ns n, since the Ns ninterface nodes are included both in Ni nandNe n. Using the piecewise linear nodal base functions Ψ q(see Appendix D), we define the finite dimensional approximations of functional spaces W1,2 01(G)3,W1,2 01(G),W1,2 02(A), and W1,2(A), W1,2,h 01(G)3=  f(r)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsinglef(r) =Ni n/summationdisplay q=Ncn+1(fx,qex+fy,qey+fz,qez) Ψq(r)  ,(3.51) W1,2,h 01(G) =  f(r)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsinglef(r) =Ni n/summationdisplay q=Ncn+1fqΨq(r)  , (3.52) W1,2,h 02(A) =  f(r)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsinglef(r) =Ni n+Ne n−Ns n−No n /summationdisplay q=Nin−Nsn+1fqΨq(r)  , (3.53) W1,2,h(A) =  f(r)/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsinglef(r) =Ni n+Ne n−Ns n/summationdisplay q=Nin−Nsn+1fqΨq(r)  , (3.54) and expand the potentialsiA,iΦ,iU, and the boundary conditioniˆUat thei-th time step into series, iA(r) =Ni n/summationdisplay q=Ncn+1/parenleftbigiAx,qex+iAy,qey+iAz,qez/parenrightbig Ψq(r), (3.55) iΦ(r) =Ni n/summationdisplay q=Ncn+1iΦqΨq(r), (3.56) 44 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH iU(r)−iˆU(r) =Ni n+Ne n−Ns n−No n /summationdisplay q=Nin−Nsn+1iUqΨq(r), (3.57) iˆU(r) =Ni n+Ne n−Ns n/summationdisplay q=Nin+Nen−Nsn−Non+1iˆUqΨq(r). (3.58) Note that the base functions corresponding to nodes placed on boundaries ∂G1and∂A2 are left out from definitions of spaces W1,2,h 01(G)kandW1,2 02(A), respectively. We arrange the coefficients of electromagnetic potentials at the i-th time step,iAx,q, iAy,q,iAz,q,iΦq, andiUq, into vector ix=/bracketleftBig/parenleftbigiAx,q,iAy,q,iAz,q,iΦq/parenrightbigNi n q=Ncn+1,/parenleftbigiUq/parenrightbigNi n+Ne n−Ns n−No n q=Nin−Nsn+1/bracketrightBig , (3.59) with dimension dx= 4 (Ni n−Nc n) +Ne n−No n, and the coefficients of the scalar magnetic potential prescribed at the outer boundary,iˆUq, into vector iy=/parenleftBig iˆUq/parenrightBigNi n+Ne n−Ns n q=Nin+Nen−Nsn−Non+1, (3.60) with dimension dy=No n. By substituting (3.55–3.58) into the integral equations (3.48– 3.50), and taking successively the base functions exΨp,eyΨp,ezΨp, Ψpforp=Nc n+ 1,...,Ni nin place of δA,δΦ, and Ψ pforp=Ni n−Ns n+ 1,...,Ni n+Ne n−Ns n−No nin place ofδU, we obtain a linear system A ·i+1x=B ·ix+C ·(iy+i+1y). (3.61) Matrices A,B, and Chave dimensions dx×dx,dx×dx, anddx×dy, respectively. Kˇ r´ ıˇ zek & Neittaanm¨ aki (1990) introduce two technique which can be used to construct the finite element matrices, “node by node” and “element by element.” We employ the later one where the matrices A,B, andCare assembled from the submatrices associated to individual tetrahedra. LetTκbe the tetrahedra of the interior ( κ= 1,...,Ni t) and exterior ( κ=Ni t+ 1,...,Ni t+Ne t) mesh. The mantle-atmosphere interface ∂G2is approximated by trian- gular facets Fη,η= 1,...,Ns twith outer normals nηoriented outward from the Earth’s centre. Let the conductivity σbe a piecewise constant function, σ(r) =σκ= const in Tκ, κ= 1,...,Ni t. (3.62) Then we construct the matrices A,B,Cas sums of contributions of individual interior tetrahedra, surface facets, and exterior tetrahedra, A=Ni t/summationdisplay κ=1(Pκ+Qκ) +Ns t/summationdisplay η=1Rη+Ni t+Ne t/summationdisplay κ=Ni t+1Sκ, (3.63) 3.1.A−Φ,UFORMULATION 45 B=Ni t/summationdisplay κ=1(Pκ− Q κ)−Ns t/summationdisplay η=1Rη−Ni t+Ne t/summationdisplay κ=Ni t+1Sκ, (3.64) C=Ni t+Ne t/summationdisplay κ=Ni t+1Tκ. (3.65) Matrices PκandQκcorresponding to the volume integrals in equations (3.48–3.49) consist respectively of 4 ×4 and 3 ×3 blocks, 4 (q−Nc n)−3 · · · · · · 4 (q−Nc n) µ0σκ ∆t Ipq κ 0 0 ex·Jpq κ 0Ipq κ 0 ey·Jpq κ 0 0 Ipq κ ez·Jpq κ ex·Jqp κey·Jqp κez·Jqp κKpq κ 4 (p−Nc n)−3 ... ... 4 (p−Nc n),(3.66) 4 (q−Nc n)−3 · · · 4 (q−Nc n)−1 1 2 Kpq κ 0 0 0Kpq κ 0 0 0 Kpq κ 4 (p−Nc n)−3 ... 4 (p−Nc n)−1,(3.67) for interior nodes Nc n< p,q ≤Ni n,rp,rq∈Tκ. Volume integrals Ipq κ,Jpq κ, andKpq κare defined in Appendix D. Matrices Rηrepresent the contribution of the surface integrals over the mantle- atmosphere interface. They include 3 ×1 and 1 ×3 blocks, 4 (Ni n−Nc n)+q−(Ni n−Ns n) −1 2 ex·Lpq η ey·Lpq η ez·Lpq η 4 (p−Nc n)−3 ... 4 (p−Nc n)−1,(3.68) 4 (q−Nc n)−3 · · · 4 (q−Nc n)−1 1 2/parenleftBig ex·Lpq ηey·Lpq ηez·Lpq η/parenrightBig 4 (Ni n−Nc n)+p−(Ni n−Ns n),(3.69) for surface nodes Ni n−Ns n<p,q ≤Ni n, and rp,rq∈Fη. See Appendix D for definition of surface integrals Lpq η. Note the opposite signs due to the orientation of outer normals on∂G2and∂A1. 46 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH Matrices Sκinclude the terms corresponding to the volume integral (3.50), 4 (Ni n−Nc n)+q−(Ni n−Ns n)/parenleftbig −1 2Kpq κ/parenrightbig 4 (Ni n−Nc n)+p−(Ni n−Ns n),(3.70) for exterior mesh nodes Ni n−Ns n<p,q ≤Ni n+Ne n−Ns n−No n, and rp,rq∈Tκ. Similarly, matrices Tκrepresenting the Dirichlet boundary value consist of terms, q−(Ni n+Ne n−Ns n−No n)/parenleftbig −1 2Kpq κ/parenrightbig 4 (Ni n−Nc n)+p−(Ni n−Ns n),(3.71) forNi n−Ns n<p≤Ni n+Ne n−Ns n−No n,Ni n+Ne n−Ns n−No n<q≤Ni n+Ne n−Ns n, and rp,rq∈Tκ. Since Pκ,Qκ,Sκare symmetric matrices, and sinceNs t/summationtext η=1Rηis symmetric due to equa- tion (D 14), the resulting Galerkin matrices AandBare also symmetric. Their lower triangular parts can be stored effectively using a sparse storage scheme (Press et al., 1996). Matrix Ais indefinite. While the first 4 ( Ni n−Nc n) diagonal terms corresponding to the interior mesh nodes are positive, the next Ne n−No ndiagonal terms corresponding to the exterior mesh nodes are negative. Various numerical methods can be applied to solve the linear system (3.61). Using the complete LU decomposition of A(IMSL, 1994) is the fastest way to solve (3.61) repeat- edly. However, the memory and computational requirements of the initial factorization ofAlimit its use only to problems with small spatial resolution. One of the methods designed for solution of large sparse symmetric linear systems is the preconditioned conjugate gradient (PCG) method (Press et al., 1992). Although the convergence is assured only for positive definite matrices, it can be used as a general iterative procedure even for indefinite symmetric systems. The incomplete LU decompo- sition, A ≈ L · U (3.72) which neglects additional fill-ins during the factorization of Ais computed initially and is used as a preconditioner for the iterative process. At each time level ti+1, the iterations are started from the solution obtained in the previous time step,ix, and the right-hand side bof (3.61) is constructed, i+1x0=ix, (3.73) b=B ·ix+C ·(iy+i+1y). (3.74) Although the actual difference betweeni+1xandixdepends on the length of time step, and on the corresponding change of the boundary condition, the initial guess (3.73) reduces the number of PCG iterations compared with the usual zero starting value. 3.2. VALIDATION 47 Then, using the auxiliary vectors rk,pk,zk, and coefficients α,βwith initial values, r0=b− A ·i+1x0, (3.75) L · U · z0=r0, (3.76) p0=z0, (3.77) we repeatedly compute for k= 1,2,..., αk=zk−1·rk−1 pk−1· A ·pk−1, (3.78) i+1xk=i+1xk−1+αkA ·pk−1, (3.79) rk=rk−1−αkA ·pk−1, (3.80) L · U · zk=rk, (3.81) βk=zk·rk zk−1·rk−1, (3.82) pk=zk+βkpk−1, (3.83) until the norm of the residuum vector rk, scaled by the norm of the right-hand side, ε=/parenleftbiggrk·rk b·b/parenrightbigg1 2 , (3.84) is smaller than a given tolerance ε0. We use the value ε0= 10−6in following computations. During each iteration step we compute one matrix multiplication A ·pk−1, and solve the factorized linear system (3.81) by backsubstitution. The number of iterations is variable, and generally increases for complicated 3-D conductivity models. We have also employed another preconditioned iterative method used by Everett & Schultz (1996). However it required about twice more iterations than the PCG method during test runs with various conductivity models. 3.2 Validation 3.2.1 Axially symmetric conductivity model Firstly we validate the 3-D finite element method against the 2-D time-domain spectral- finite element approach by Martinec et al. (2002) which computes the response of axially symmetric conductivity models to axially symmetric transient external excitation sig- nals. The 2-D solution does not represent an exhaustive test of 3-D methods, since it does not include the galvanic term µ0σgrad∂Φ ∂tin the equation of EM induction (3.21), 48 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH Figure 3.2: The longitudinal component of magnetic vector potential, Aϕ, in nA for the axially sym- metric model excited by Load 2 shown for two time levels. Results obtained by the 3-D finite element computation are plotted on the left hemisphere. The right hemisphere corresponds to the 2-D spectral- finite element time-domain approach, with highly conductive core ( σ= 1014S/m). and can be expressed using only the longitudinal (toroidal) component of the vector mag- netic potential Aϕ(Everett & Schultz, 1996; Martinec et al., 2002). On the other hand, the potentials Aϕobtained by both methods are directly comparable, and subsequent numerical differentiation to compute Bis not necessary to perform. In the presented example we prescribe a homogeneous mantle with conductivity 1S/m. In the northern hemisphere we overlay the mantle by a 700 km thick layer with conduc- tivity 0.1S/m. Since the 2-D spectral-finite element method does not assume an infinitely conductive core, we approximate it by a sphere with conductivity of 1014S/m. The model is excited by the analytical model of a geomagnetic storm introduced as “Load 2” by equa- tion (2.111) in Chapter 2, with amplitude A= 0.003nT/sand relaxation time 1 /α= 48 h. The external excitation field is prescribed directly by the means of the dipolar spherical harmonic coefficient in the 2-D method. In the 3-D finite element method, the Dirichlet boundary condition prescribed at nodes rq∈∂A2is constructed by the substitution of (2.111) into (2.11), iˆUq=Ate−α tzq. (3.85) Figure 3.2 compares two snapshots of Aϕobtained by both methods at two time levels. The 3-D finite element solution was computed on a regular tetrahedral mesh 3.2. VALIDATION 49 Ar Aϑ Aϕ -8-6-4-2024(µA) 0200400 t (h)-8-6-4-2024(µA) 0200400 t (h)-8-6-4-2024(µA) 0200400 t (h)-30-20-10010(µA) 0200400 t (h)-30-20-10010(µA) 0200400 t (h)-30-20-10010(µA) 0200400 t (h)-200-150-100-500(mA) 0200400 t (h)-200-150-100-500(mA) 0200400 t (h)-200-150-100-500(mA) 0200400 t (h) Figure 3.3: Time evolution of the spherical components of Aat the surface point ( ϑ, ϕ) = (45◦,0) is shown for the axially symmetric model excited by Load 2. Dashed, dotted, and solid lines correspond to the 3-D finite element solutions on meshes M5,20,M5,33, and to the 2-D spectral-finite element solution, respectively. M5,20(see Table C.1). The 2-D spectral-finite element method used 100 finite elements in the radial direction, and the spherical harmonic expansion was truncated at jmax= 40. The integration time step in both methods was set to 1 h. The time-evolution of both solution at a particular point at the surface is plotted in Figure 3.3. Good agreement is obtained in both the spatial and temporal shape of Aϕ. The finite triangulation of the Earth’s surface causes that the axisymmetry of the 3-D finite element solution is disturbed which leads to small non-zero components ArandAϑ. Moreover, Figure 3.3 shows also the 3-D finite element solution computed on a tetrahedral mesh M5,33. While this mesh discretizes the spherical shell Gin the same way as M5,20, the outer boundary ofA,∂A2is positioned at aO= 10.245ainM5,33, compared to aO= 4.068ain the case ofM5,20. Note that the increase of aOis not reflected in Aϕ. This validates the choice of aO≈4aused in the presented runs. 3.2.2 Nested-sphere conductivity model The next example validates the 3-D finite element solver against the semi-analytical solution for a nested sphere conductivity model introduced in Section 2.2.1. The off-axis position of the spherical inclusion yields a fully 3-D solution including the galvanic effect. The nested spheres conductivity Models A and B that were used in the validation examples in Chapter 2 contain a spherical inclusion of a comparable size to the Earth’s core. Such inclusion is too large to fit into the spherical shell Gthat approximates the Earth’s mantle in the 3-D finite element solver. Therefore we introduce another nested spheres configuration, Model C, with a smaller spherical inclusion positioned in 50 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH Figure 3.4: The high-resolution, 3-D finite element solution for nested spheres conductivity Model C excited by Load 2. Spherical components of vector magnetic potential Aand scalar electrical potential Φ are shown at t= 48 h after the onset of the storm model. The scales of Aand Φ are in nA and nV/s, respectively. Boundaries of the core and spherical inclusion are outlined. the mantle. The parameters of the model, as denoted in Figure 2.3, are a= 6371 km, b= 1200 km, d= 5000 km, ϑd= 40◦,ϕd= 30◦,σ1= 1S/m, andσ2= 10S/m. The semi-analytical method does not allow to include another disjunctive sphere that would represent a highly conductive core. Nevertheless, the choice of highly conductive host sphere prevents the EM wave to penetrate down to the core and makes both methods 3.2. VALIDATION 51 Br Bϑ Bϕ -30-20-10010(nT) 400050006000-30-20-10010(nT) 400050006000-30-20-10010(nT) 400050006000-30-20-10010(nT) 400050006000-30-20-10010(nT) 400050006000020406080100120 400050006000020406080100120 400050006000020406080100120 400050006000020406080100120 400050006000020406080100120 4000500060000246 4000500060000246 4000500060000246 4000500060000246 4000500060000246 400050006000 -30-20-10010(nT) 400050006000-30-20-10010(nT) 400050006000-30-20-10010(nT) 400050006000-30-20-10010(nT) 400050006000-30-20-10010(nT) 400050006000020406080100120 400050006000020406080100120 400050006000020406080100120 400050006000020406080100120 400050006000020406080100120 4000500060000246 4000500060000246 4000500060000246 4000500060000246 4000500060000246 400050006000 -30-20-10010(nT) 400050006000 r (km)-30-20-10010(nT) 400050006000 r (km)-30-20-10010(nT) 400050006000 r (km)-30-20-10010(nT) 400050006000 r (km)-30-20-10010(nT) 400050006000 r (km)020406080100120 400050006000 r (km)020406080100120 400050006000 r (km)020406080100120 400050006000 r (km)020406080100120 400050006000 r (km)020406080100120 400050006000 r (km)0246 400050006000 r (km)0246 400050006000 r (km)0246 400050006000 r (km)0246 400050006000 r (km)0246 400050006000 r (km) Figure 3.5: 3-D finite element solution for the nested spheres conductivity Model C excited by Load 2. The spherical components of Bobtained by MLS fit are plotted along radius at fixed colatitude ϑ= 30◦and longitude ϕ= 45◦. Open circles and dots correspond to low-resolution mesh Mr 4,10and high-resolution mesh Mr 5,20, respectively. For comparison, the semi-analytical solution and the solutions obtained by the spectral-finite element approach introduced in Chapter 2, are shown by solid and dashed lines, respectively. Top, middle, and bottom plates correspond to time levels 24, 48, and 72 h after the onset of the analytical storm model. 52 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH Br Bϑ Bϕ -30-20-100(nT) 0200400 t (h)-30-20-100(nT) 0200400 t (h)-30-20-100(nT) 0200400 t (h)-30-20-100(nT) 0200400 t (h)020406080100120140 0200400 t (h)020406080100120140 0200400 t (h)020406080100120140 0200400 t (h)020406080100120140 0200400 t (h)024 0200400 t (h)024 0200400 t (h)024 0200400 t (h)024 0200400 t (h) Figure 3.6: Same result as in Figure 3.5, but now we show the time evolution of Bat the surface point (ϑ, ϕ) = (30◦,45◦). comparable. The conductivity Model C is again excited by Load 2, the analytical approximation of a geomagnetic storm. The semi-analytical solution was computed for discrete frequencies from the interval 10−9–10−4rad/s, and transformed to the time-domain by the inverse Fourier transformation. The truncation degree of the spherical harmonic expansion varies 0100200300400t (h) 0200400600 k -6-4-202log(ε) Figure 3.7: Rate of convergence of the preconditioned conjugate gradient solver for the high-resolution solution of the nested spheres model. Length of each horizontal line corresponds to the number of iterations of the method at given time level, while the residual norm εis shown by colors. 3.2. VALIDATION 53 from 12 to 24, according to the stability and accuracy criteria discussed in Section 2.2.1. The 3-D finite element solution is found for two tetrahedral meshes, a low-resolution meshMr 4,10and a high-resolution mesh Mr 5,20(see Table C.1). Both meshes are locally refined inside the spherical inclusion and in its neighbourhood. The integration time step is 1 h. The magnetic induction vector Bis obtained by a numerical differentiation of potential Ausing the Moving Least Squares (MLS) method described in Appendix E. Figure 3.4 shows the snapshot of the solution at t= 48 h in terms of the potentials A, Φ. The anomalies of Ar,Aϑ, and Φ are observed near the spherical inclusion, and to a lesser extent near the Earth’s surface, where the axial symmetry is also disturbed. Figures 3.5 and 3.6 compare the spatial and temporal behaviour of the 3-D finite element solutions with the semi-analytical reference solution. We also display the response of Model C computed by the 3-D spectral-finite element method introduced in Chapter 2, at a high spatial resolution ( jmax= 40,P= 100, ∆t= 1 h). The results indicate that the low-resolution 3-D finite element solution is accurate in terms of BrandBϑcomponents, but yields large discrepancies in the Bϕcomponent. The accuracy of the high-resolution solution is comparable to that reached by the spectral-finite element approach. Finally, Figure 3.7 shows the rate of convergence of the PCG method for the high- resolution 3-D finite element computation. As expected, the largest number of iterations per time step corresponds to the areas of steepest increase or descent of the inducing load. 54 3. TIME-DOMAIN, 3-D FINITE ELEMENT APPROACH 55 Chapter 4 The transient Dst induction signal at satellite altitudes 4.1 Introduction In this chapter, we consider a realistic crust and upper-mantle conductivity model based on laboratory measurements and seismic tomography. This model is excited by a tran- sient Dst signal taken from the period 1979–1980. Using the time-domain, spectral-finite element approach introduced in Chapter 2 we will study the sensitivity of the EM response observed at satellite altitudes to lateral conductivity variations in the upper mantle which are overlain by the near-surface conductivity inhomogeneities characterizing the ocean, continent and sea-shelf distribution. 4.2 Conductivity model We construct a 3-D conductivity model based on three data sources. The contrasts between the highly conductive oceans water, resistive continents, and intermediate shelves and marine sediments are treated using the surface conductance map constructed by Everett et al. (2002). The conductance is converted to conductivity using a layer with constant thickness of 50 km (see Figure 4.1). The conductivity in the upper and lower mantle is based on the 1-D conductivity model derived by Xu et al. (2000) from laboratory measurements, combined with the seismic tomography model SKS12WM13 by Liu & Dziewonski (1994). Although the shear-wave velocity and electrical conductivity have different sensitivity to various physi- cal and chemical properties of the mantle, the positive correlation between the conductor depths and velocity perturbations obtained by Tarits (1994) supports the basic corre- spondence between the fast, cold, and resistive mantle contrasting to the slow, hot, and conductive areas. Approximating the relation between electrical conductivity and shear-wave velocity variationδVby an exponential function, the conductivity σat depthh, colatitude ϑ, and longitudeϕis σ(h,ϑ,ϕ ) =σ0(h)f(h)−δV(h,ϑ,ϕ ) δVM(h)−δVm(h), (4.1) 564. THE TRANSIENT DST INDUCTION SIGNAL AT SATELLITE ALTITUDES 0 200 400 600 800 1000h (km) -4-3-3-2-101 log(σ0 (S/m))0 200 400 600 800 1000h (km) -4-3-3-2-101 log(σ0 (S/m))0 200 400 600 800 1000h (km) -4-3-3-2-101 log(σ0 (S/m)) -3.0-2.5-2.0-1.5-1.0-0.5 -3.0-2.5-2.0-1.5 0.00.5400 km 800 km0−50 km -3.0-2.5-2.0-1.5 -1.5-1.0-0.50.00.5 0.00.5200 km 600 km 1000 km Figure 4.1: 3-D conductivity model in the crust, upper and mid-mantle. The solid line in the left panel shows the 1-D profile based on laboratory measurements (Xu et al., 2000). Lateral variations are confined to the interval shown by dotted lines. The first color map on the right shows the conductivity obtained from the surface conductance map for a uniform depth of 50 km. The following cross-sections sample the tomography-derived conductivity at depths of 200, 400, 600, 800, and 1000 km. All scale bars are in units of log( σ(S/m)). whereδVM(h) andδVm(h) are respectively the maximum and minimum lateral seismic- velocity variations at depth htaken from SKS12WM13, and σ0(h) is the 1-D conductivity model by Xu et al. (2000). The base f(h) which is the ratio of the maximum to minimum conductivity in each layer, is set to 100 in the upper mantle and to 10 in the lower mantle (Shankland et al., 1993). Finally, a homogeneous core with conductivity 10000S/mis assumed. In order to evaluate the compatibility of the 3-D conductivity model with surface 4.2. CONDUCTIVITY MODEL 57 -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)ABN 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)ALM 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)API -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)AML 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)HON -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)CLH 584. THE TRANSIENT DST INDUCTION SIGNAL AT SATELLITE ALTITUDES 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)MBO -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)KNY0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)KAK -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)IRT 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)PAG -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)MMB 4.2. CONDUCTIVITY MODEL 59 -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)PIL 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)PRU 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)SJG -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)RSV 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)SSO -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)SOD 604. THE TRANSIENT DST INDUCTION SIGNAL AT SATELLITE ALTITUDES -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)TOL 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)TUC 0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)0.010.010.02 0.050.10.2 f (cpd)WIT -5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)-5000500100015002000C (km) 0.010.010.02 0.050.10.2 f (cpd)WAT Figure 4.2: Comparison of the Schmucker C-responses generated by the 3-D conductivity model with realistic responses computed by Schultz & Larsen (1987) for 22 ground stations. Response of the 1- D conductivity model σ0(h) which is identical at all observatories is plotted by dashed lines, while the responses of the 3-D conductivity model are shown by solid lines. The observed observatory responses are confined within the one standard deviation error shown by color bands. Red, and blue colors correspond to Real( C) and −Imag( C), respectively. geomagnetic data, the Schmucker C-responses, C=atanϑBr 2Bϑ, (4.2) are computed in the frequency range 0 .01 to 0.2 cpd using the spectral-finite element code by Martinec (1999a). Figure 4.2 shows them compared with local 1-D responses from the 4.2. CONDUCTIVITY MODEL 61 Station Geomagnetic Geographic χ2misfit latitude latitude longitude 1-D 3-D ABN 53.65 51 .19 359 .61 310.33 570.87 ALM 40.23 36 .85 357 .54 78.85 62.03 AML −47.25 −43.15 172 .72 696.78 1409.18 API −15.69 −13.81 188 .23 117.40 261.23 CLH 49.85 38 .73 283 .16 378.98 308.43 HON 21.41 21 .32 201 .94 150.09 520.40 IRT 41.30 52 .46 104 .04 97.22 99.59 KAK 26.44 36 .23 140 .19 102.96 102.52 KNY 20.93 31 .42 130 .88 71.40 55.42 MBO 20.80 14 .39 343 .04 17.92 67.60 MMB 34.44 43 .91 144 .19 77.27 830.06 PAG 40.60 42 .52 24 .18 385.67 112.98 PIL −20.54 −31.67 296 .12 131.51 19.03 PRU 49.61 49 .99 14 .55 136.75 120.08 RSV 55.56 55 .84 12 .46 88.64 62.69 SJG 29.54 18 .38 293 .88 434.86 153.28 SOD 63.62 67 .37 26 .65 223.28 200.20 SSO 23.45 33 .58 135 .94 191.18 167.04 TOL 43.48 39 .88 355 .95 336.32 315.25 TUC 40.47 32 .25 249 .17 151.94 345.32 WAT −41.44 −30.32 115 .88 404.71 193.02 WIT 53.82 52 .81 6 .67 111.01 132.07 Table 4.1: Misfit of the C-responses generated by the 1-D conductivity model σ0(h) (Xu et al., 2000), and by the 3-D conductivity model (4.1), respectively, with respect to the observatory values from the Schultz & Larsen (1987) dataset. Blue and red colors are respectively used where the introduction of lateral variations improves or deteriorates the fit of the data. dataset by Schultz & Larsen (1987). The agreement between the synthetic and observed data at individual stations is highly variable with the χ2misfit ranging from 19 .03 up to 1409.18 (see Table 4.1). Addition of the lateral conductivity variations according to (4.1) to the 1-D profile σ0(h) improved the fit to the observed C-responses at 13 stations. At the remaining 9 stations the misfit obtained for the 3-D model is larger than for the 1-D model. These results must be interpreted cautiously, since the 3-D conductivity model does not comply with the assumption of local 1-D structure used in the processing of local responses. Nevertheless, stacking over the computed and observed mid-latitude C- responses (see Figure 4.3), shows that the model has a degree of lateral heterogeneity appropriate to the observed responses, although the particular conductivity structures in 624. THE TRANSIENT DST INDUCTION SIGNAL AT SATELLITE ALTITUDES -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd)050010001500Real C (km) -1000-5000Imag C (km) 0.010.010.02 0.050.10.2 frequency (cpd) Figure 4.3: The Schmucker C-responses generated by the 3-D conductivity model. At middle geo- magnetic latitudes, 25◦≤ |λd| ≤65◦, these are limited to the areas shown by grey shading. Local 1-D responses computed by Schultz & Larsen (1987) for 22 ground stations are stacked over using solid lines. the model may differ from those required by the surface data. 4.3 Exciting field We excite the conductivity model specified in the previous section by a transient, axisym- metric ring current. This is a first-order approximation of the spatial distribution of Dst -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h)1983 -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h)1990 -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h)2020 -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h)2229 -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h)2253 -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h)2441 -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h)2457 -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h)2480 -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h)2618 -1000100200G10(e) (nT) 2000220024002600 t (h)-1000100200G10(e) (nT) 2000220024002600 t (h) Figure 4.4: The coefficient G(e) 10of the external excitation field. Time t= 0 corresponds to Nov. 5, 1979, 0:00 UT, the displayed interval spans over three successive storms occurring between Jan. 24 and Feb. 27, 1980. Stars mark the epochs shown in Figures 4.5–4.7. 4.4. RESULTS 63 currents (Daglis & Kozyra, 2002). Assuming a P10geometry for the exciting field in the geomagnetic coordinate system, the dipolar coefficient of the external magnetic potential G(e) 10can be derived from the Dst index (downloaded from the WDC-C2 Kyoto Dst index service, http://swdcdb.kugi.kyoto-u.ac.jp/dstdir/ ), by a simple scaling formula, G(e) 10(t) =−/radicalbigg 4π 3Dst(t) 1 +Q1, (4.3) whereQ1= 0.27 is the global estimate of the internal to external fields ratio (Langel & Estes, 1985). Figure 4.4 shows the time evolution of G(e) 10during three successive magnetic storms in January and February, 1980. At this epoch the solution has already evolved for more than 1900 h to ensure that the presented results are not biased by the choice of initial value of the magnetic field, which is an inherent feature of the time-domain method. 4.4 Results The time-domain response of the realistic 3-D conductivity model to the Dst excitation is computed using the spectral-finite element approach described in Chapter 2. The spherical harmonic expansion is truncated at jmax= 40, yielding an angular resolution of about 500 km at the Earth’s surface. The radial discretization assumes layer thickness of 25 km throughout the mantle, and 200 km in the core, and the time-step is set to 1 h. Figures (4.5–4.7) show the snapshots of the Br,Bϑ, andBϕcomputed at 400 km altitude which is typical of low-orbit satellites. Note that the spherical components are related to the geomagnetic coordinate system. In the case of an azimuthally symmet- ric external excitation and a 1-D conductivity model, the longitudinal component Bϕ vanishes. Therefore, Bϕshown in Figure (4.7) can be completely assigned to the effect of lateral conductivity variations. On the other hand, the BrandBϑcomponents have predominantly dipolar shape with amplitudes reaching 50 nT, and 130 nT, respectively, at the peak of the most intense storm. About 50% of Br, and 30% of Bϑat this altitude is of internal origin. In order to make the effect of lateral variations on Br, andBϑvisible, both the internal, and external dipolar terms are subtracted from the solution in Figures 4.5 and 4.6. The largest magnetic anomalies, up to 1.2 nT, 1.8 nT, and 2.1 nT in Bϕ, and in non-dipolar parts of Bϑ,Br, respectively can be assigned to the resistive areas in the transition zone below South America, Africa, and the western Pacific. On the other hand, the signal due to the large positive anomaly at 200 km depth below the eastern Pacific is weaker. The induced field shows similar spatial structure for different storms with the largest amplitude corresponding to the most severe storm. This is not surprising, since the imposed inducing field has the same dipolar spatial structure all the time. Note 644. THE TRANSIENT DST INDUCTION SIGNAL AT SATELLITE ALTITUDES -2-1012Br (nT) 1983 hrs 2229 hrs 2457 hrs 2480 hrs2253 hrs1990 hrs 2618 hrs2441 hrs2020 hrs Figure 4.5: The computed Branomalies of the induced field at satellite altitude 400 km. Only non- dipolar part of the field is shown, the terms corresponding to the external excitation, and to the induced P10field are not included. Snapshots at t= 1983, and 2441 h correspond to the onsets of the magnetic storms. At t= 1990, 2020, 2229, 2457 h the storms culminate, and t= 2253, 2480, and 2618 h sample the recovery phases. The time snapshots are marked in 4.4. that in the last snapshots in Figures (4.5–4.7), the small positive Dst index (i.e., negative G(e) 10) yields change of the sign also in the induced anomalous field. 4.5 Discussion Our simulation shows that the EM response of the realistic conductivity model to in- tense transient changes in the ring-current system is detectable by a low-orbit satel- lite magnetometer. The magnitude of the signal induced by the highly heterogeneous oceanic/continental layer is less than 10% smaller than that of the mid-mantle conduc- tivity heterogeneities. This is in agreement with the results of previous synthetic studies 4.5. DISCUSSION 65 -101Bϑ (nT) 1983 hrs 2229 hrs 2457 hrs 2480 hrs2253 hrs1990 hrs 2618 hrs2441 hrs2020 hrs Figure 4.6: The computed Bϑanomalies of the induced field at satellite altitude 400 km. See also Figure 4.5 caption. carried out in the frequency domain using surface observations (Weiss & Everett, 1998). Since the period we covered coincides with the MAGSAT mission, direct comparison of the results is advisable. Constable et al. (manuscript, http://mahi.ucsd.edu/Ste- ve/MDAT/ , 2001) isolated the Dst signal from MAGSAT data, and showed that after fitting theP10model for individual satellite passes, the RMS residuals are 5 to 6 nT for Br, andBϑcomponents, and up to 11 nT in the Bϕcomponent for the ascending passes of the spacecraft. These residuals represent, beside possible remnants of non-Dst signals, a combined effect of non- P10external sources, and lateral conductivity variations. In order to distinguish between these two components and to interpret satellite magnetic data in terms of laterally varying conductivity models, a more accurate model of the temporal andspatial variations of the external excitation field will be necessary. The combination of satellite data from several satellites and land-based measurements should allow to go beyond the simple P10excitation model. The time resolution needed to resolve rapid 664. THE TRANSIENT DST INDUCTION SIGNAL AT SATELLITE ALTITUDES -1.0-0.50.00.51.0Bϕ (nT) 1983 hrs 2229 hrs 2457 hrs2253 hrs1990 hrs 2480 hrs2020 hrs 2441 hrs 2618 hrs Figure 4.7: The computed Bϕanomalies of the induced field at satellite altitude 400 km. See also Figure 4.5 caption. changes of the external field is of the order of a typical orbital period of a satellite ( ∼1 h). 67 Chapter 5 Conclusions The time-domain approaches to the problem of electromagnetic induction in the het- erogeneous Earth’s mantle presented in Chapters 2 and 3 represent alternative methods to the traditional frequency-domain techniques. Both approaches are able to compute the response of 3-D conductivity models to external excitation with complicated spatio- temporal characteristics. The methods have been numerically validated against semi- analytical solutions for nested spheres conductivity models and numerical 2-D axially symmetric solutions. The spectral-finite element method described in Chapter 2 offers an efficient treatment of the boundary conditions at the Earth’s surface based on the analytical expression for the magnetic potential in the atmosphere. Using equations (2.13–2.15), the magnetic field computed at the Earth’s surface can be easily continued upwards to be compared with data recorded by satellite born magnetometers. Implicit treatment of the effect of 1-D spherically symmetric conductivity stabilizes the time-integration scheme, while the computationally demanding coupling due to the lateral conductivity variations is efficiently evaluated on the angular grid using the solution from the previous time step. The spectral-finite element spatial parameterization allows one to increase the radial resolution where necessary, e.g., in the lithosphere, or throughout the transition zone in the mantle. On the other hand, the angular resolution given by the truncation degree of the spherical harmonic expansion can be increased only globally. The time step has to be chosen cautiously with respect to the vigour of the time-variations of the inducing magnetospheric currents. The 3-D finite element method introduced in Chapter 3 addresses the problem of non- uniform spatial distribution of ground geomagnetic observatories. Its main advantage over the spectral-finite element method is the versatility of the spatial discretization. Using the local refinement technique, the mesh can be adapted to yield high-resolution results for regions with good data coverage. On the other hand, the 3-D finite element formulation with a prescribed model of external inducing currents requires to compute the numerical solution of the Laplace equation in the insulating atmosphere, introducing an additional computational burden. The speed of the 3-D finite element method depends on the rate of convergence of the PCG algorithm can be fasten up using a better preconditioner. This can be exploited in the inverse problem, where a complete LU decomposition can be computed for a reference 68 5. CONCLUSIONS 3-D conductivity model and used as a preconditioner for a set of forward problem runs. The presented approaches to the problem of EM induction in the Earth’s mantle are designed to be used as forward solvers in the global inversion of the satellite, or combined satellite and ground based data in the time-domain, to reveal the lateral conductivity het- erogeneities in the mid-mantle. As a preliminary step, the spectral-finite element method was applied to compute the satellite altitude response of a realistic mantle conductivity model to a Dst transient signal (see Chapter 4). The results suggest that the signal due to the mid-mantle lateral conductivity variations induced by intense geomagnetic storms should be detectable at satellite altitudes. REFERENCES 69 References Anderson, E., Z. Bai, C. Bischof, S. Blackford, J. Dongara, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, & D. Sorensen, 1999. LAPACK User’s Guide, Third edition, SIAM, Philadelphia, Pennsylvania. Banks, R.J., 1969. Geomagnetic variations and the electrical conductivity of the upper mantle, Geophys. J. R. astr. Soc., 17, 457–418;2 87. Biro, O. & K. Preis, 1989. On the use of the magnetic vector potential in the finite element analysis of three-dimensional eddy currents, IEEE Trans. Magn., 25 3145–3159. Biro, O. & K. Preis, 1990. Finite element analysis of 3-D eddy currents, IEEE Trans. Magn., 26 418–423. Daglis, I.A. & J.U. Kozyra, 2002. Outstanding issues of ring current dynamics, Journal of Atmospheric and Solar-Terestrial Physics, 64, 253–264. Duffy, T.S. & R.J. Hemley, 1995. Some like it hot: The temperature structure of the Earth, Rev. Geophys. Vol. 33 Suppl., American Geophysical Union. Everett, M.E. & A. Schultz, 1995. Geomagnetic induction in eccentrically nested spheres, Phys. Earth planet. Inter., 92, 189–198. Everett, M.E. & A. Schultz, 1996. Geomagnetic induction in a heterogeneous sphere: Azimuthally symmetric test computations and the response of an undulating 660-km discontinuity, J. geophys. Res., 101, B2, 2765–2783. Everett, M.E., 1997. A three-dimensional spherical mesh generator, Geophys. J. Int., 130, 193–200. Everett, M.E., S. Constable, & C.G. Constable. Effect of near-surface conductance on global satellite induction responses, submitted to Geophys. J. Int. Fainberg, E.B., A.V. Kuvshinov, & B.Sh. Singer, 1990. Electromagnetic induction in a spherical earth with non-uniform oceans and continents in electric contact with the underlying medium — I. Theory, method and example, Geophys. J. Int., 102, 273–281. Gossier, A., 2001. Moving least-squares: A numerical differentiation method for irreg- ularly spaced calculation points, Sandia report SAND2001-1669, Albuquerque, New Mexico and Livermore, California. 70 REFERENCES Hamano, Y., 2002. A new time-domain approach for the electromagnetic induction prob- lem in a three-dimensional heterogeneous earth, Geophys. J. Int., 150, 753–769. IMSL Fortran subroutines for mathematical applications, 1994. Math library, Volumes 1 & 2, version 3.0, Visual Numerics, Inc., Houston, Texas. Kˇ r´ ıˇ zek, M. & P. Neittaanm¨ aki, 1990. Finite Element Approximation of Variational Prob- lems and Applications, Longman Scientific and Technical, Harlow, UK and Wiley, New York. Kuvshinov, A.V. & O.V. Pankratov, 1994. Electromagnetic induction in a spherical Earth with inhomogeneous conducting mantle: thin sheet forward problem, in: Proc. 12th workshop on electromagnetic induction in the Earth, p. 53, Universit´ e de Bretagne Occidentale, Brest, France. Lahiri, B. & A. Price, 1939. Electromagnetic induction in non-uniform conductors, and the determination of the conductivity of the Earth from terrestrial magnetic variations, Phil. Trans. R. Soc. London, 237, 509–540. Langel, R.A. & R.H. Estes, 1985. Large-scale, near-field magnetic fields from external sources and the corresponding induced internal field, Geophys. J. Int., 90, 2487–2494. Langel, R.A., 1987. The main field, in: Geomagnetism, Volume 1, ed. Jacobs J.A., Aca- demic Press, London. Langel, R.A., T.J. Sabaka, R.T. Baldwin, & J.A. Conrad, 1996. The near-Earth magnetic field from magnetospheric and quiet-day ionospheric sources and how it is modeled, Phys. Earth planet. Inter., 98, 235–267. Liu, X.-F. & A.M. Dziewonski, 1994. Improved resolution of the lower-most mantle shear wave velocity structure obtained using SKS-S data, EOS, Trans. Am. geophys. Un., 75, Spring Meet. Suppl., 232–234. Liu, A. & B. Joe, 1996. Quality local refinement of tetrahedral meshes based on 8- subtetrahedron subdivision, Math. Comp., 215 1183-1200. Martinec, Z., 1989. Program to calculate the spectral harmonic expansion coefficients of the two scalar fields product, Computer Physics Communications, 54, 177–182. Martinec, Z., 1997. Spectral-finite element approach to two-dimensional electromagnetic induction in a spherical Earth, Geophys. J. Int., 130, 583–594. Martinec, Z., 1998. Geomagnetic induction in multiple eccentrically nested spheres, Geo- phys. J. Int., 132, 96–110. REFERENCES 71 Martinec, Z., 1999. Spectral-finite element approach to three-dimensional electromagnetic induction in a spherical Earth, Geophys. J. Int., 136, 229–250. Martinec, Z., 1999. Spectral, initial value approach for viscoelastic relaxation of a spher- ical earth with a three-dimensional viscosity — I. Theory, Geophys. J. Int., 137, 469– 488. Martinec, Z., M.E. Everett & J. Vel´ ımsk´ y, 2002. Time-domain, spectral-finite element approach to transient two-dimensional geomagnetic induction in a spherical heterogoe- neous earth, submitted to Geophys. J. Int. N´ edelec, J.C., 1980. Mixed finite elements in R3,Numer. Math. 35, 315–341. Olsen, N., 1998. The electrical conductivity of the mantle beneath Europe derived from C-responses from 3 to 720 hr, Geophys. J. Int., 133, 298–308. Parkinson, W.D. & V.R.S. Hutton, 1989. The electrical conductivity of the Earth, in: Geomagnetism, Volume 3, ed. Jacobs J.A., Academic Press, London. Petersons, H.F. & R.S. Anderssen, 1990. On the spherical symmetry of the electrical conductivity of the Earth’s mantle, J. Geomag. Geoelectr., 42, 1309–1324. Press, W.H., S.A. Teukolsky, W.T. Vetterling & B.P. Flannery, 1992. Numerical recipes in Fortran. The art of scientific computing, Cambridge University Press, Cambridge. Press, W.H., S.A. Teukolsky, W.T. Vetterling & B.P. Flannery, 1996. Numerical recipes in Fortran 90. The art of parallel scientific computing, Cambridge University Press, Cambridge. Price, A.T., 1930. Electromagnetic induction in a conducting sphere, Proc. London Math. Soc., 51, 217. Rikitake, T., 1950. Electromagnetic induction within the Earth and its relation to the electrical state of the Earth’s interior, Bull. Earthq. Res. Inst., Tokyo Univ., 28, 45–100. Roberts, R.G., 1984. The long-period electromagnetic response of the Earth, Geophys. J. R. astr. Soc., 78, 547–572. Schultz, A., 1990. On the vertical gradient and associated heterogeneity in mantle elec- trical conductivity. Phys. Earth planet. Inter., 64, 68–86. Schultz, A. & J. Larsen, 1987. On the electrical conductivity of the mid mantle: I. Cal- culation of equivalent scalar magnetotelluric response functions, Geophys. J. R. astr. Soc., 88, 733–761. 72 REFERENCES Schultz, A. & J. Larsen, 1990. On the electrical conductivity of the mid mantle: II. Delin- eation of heterogeneity by application of extremal inverse solutions, Geophys. J. Int., 101, 565–580. Semenov, V.Yu., 1998. Regional conductivity structures of the Earth’s mantle, Publi- cations of the Institute of Geophysics, Polish Academy of Sciences, Vol. C-65 (302), Warsaw. Shankland T.J., J. Peyronneau, & J.P. Poirier, 1993. Electrical conductivity of the Earth’s lower mantle, Nature, 366, 453–455. Tarits, P., 1994. Electromagnetic studies of global geodynamic processes, Surveys in Geophysics, 15, 209–238. Uyeshima, M. & A. Schultz, 2000. Geoelectromagnetic induction in a heterogeneous sphere: a new three-dimensional forward solver using a conservative staggered-grid finite difference method, Geophys. J. Int., 140, 636–650. Varshalovich, D.A., A.N. Moskalev & V.K. Khersonskii, 1989. Quantum Theory of An- gular Momentum, World Scientific, Singapore. Vel´ ımsk´ y, J., M.E. Everett & Z. Marinec, 2003. The transient Dst electromagnetic induc- tion signal at satellite altitudes for a realistic 3-D electrical conductivity in the crust and mantle, Geophys. Res. Lett., in press. Weiss, C.J. & M.E. Everett, 1998. Geomagnetic induction in a heterogeneous sphere: Fully three-dimensional test computations and the response of a realistic distribution of oceans and continents, Geophys. J. Int., 135, 650–662. Xu, Y., B.T. Poe, T.J. Shankland, & D.C. Rubie, 1998. Electrical conductivity of olivine, wadsleyite, and ringwoodite under upper-mantle conditions. Science, 280, 1415–1418. Xu, Y., T.J. Shankland, & B.T. Poe, 2000. Laboratory-based electrical conductivity in the Earth’s mantle. J. geophys. Res., 105(B12), 27,865–27,875. Yoshimura, R. & N. Oshiman, 2002. Edge-based finite element approach to the simulation of geolectromagnetic induction in a 3-D sphere, Geophys. Res. Lett., 29, No. 3, 9-1–9-4. Zhang, T.S. & A. Schultz, 1992. A 3-D perturbation solution for the EM induction problem in a spherical earth — the forward problem, Geophys. J. Int., 111, 319–334. 73 Appendix A Spherical harmonic functions This appendix summarizes the definitions of the scalar and vector spherical harmonic functions which are used in Chapter 2, including normalization, derivatives, and separa- tion of radial and angular components. We expand scalar functions into series of fully normalized scalar spherical harmonic functionsYjm(Varshalovich et al., 1989), Yjm(Ω) =/radicalBigg 2j+ 1 4π(j−m)! (j+m)!Pm j(cosϑ)eimϕ=Pjm(cosϑ)eimϕ, (A 1) wherePjm(cosϑ) is the normalized form of associated Legendre polynomials Pm j(cosϑ). The spherical harmonic functions are orthonormal on a unit sphere, /integraldisplay ΩYjm(Ω)·Yj/primem/prime(Ω)dΩ =δjj/primeδmm/prime, and are eigenfunctions of the angular part of the the Laplace operator, /bracketleftbigg1 sinϑ∂ ∂ϑ/parenleftbigg sinϑ∂ ∂ϑ/parenrightbigg +1 sin2ϑ∂2 ∂ϕ2/bracketrightbigg Yjm(Ω) = −j(j+ 1)Yjm(Ω). (A 2) Similarly, vector functions can be expanded into series of vector spherical harmonics. From several different sets of orthogonal vector spherical harmonic functions, we prefer the one used by Martinec (1999b), /braceleftBig S(λ) jm(Ω)|j= 0,...,∞, m=−j,...,j, λ =−1,0,1/bracerightBig , (A 3) where S(0) jm(Ω) = er×gradΩYjm(Ω) =/bracketleftbigg∂Yjm(Ω) ∂ϑeϕ−1 sinϑ∂Yjm(Ω) ∂ϕeϑ/bracketrightbigg ,(A 4) S(−1) jm(Ω) =Yjm(Ω)er, (A 5) S(1) jm(Ω) = gradΩYjm(Ω) =/bracketleftbigg∂Yjm(Ω) ∂ϑeϑ+1 sinϑ∂Yjm(Ω) ∂ϕeϕ/bracketrightbigg . (A 6) 74 A. SPHERICAL HARMONIC FUNCTIONS The angular part of the gradient operator is gradΩ=eϑ∂ ∂ϑ+eϕ1 sinϑ∂ ∂ϕ, (A 7) ander,eϑ, and eϕdenote unit vectors corresponding to the radial coordinate r, colati- tudeϑ, and longitude ϕ, respectively. Functions S(0) jm(Ω) are toroidal, functions S(−1) jm(Ω) andS(1) jm(Ω) represent the radial and angular component of the spheroidal part of the expanded field, respectively. Although S(λ) jm(Ω) are not eigenfunction of the Laplace op- erator, separation of the radial component generally simplifies the evaluation of coupling in problems with laterally varying parameters. Functions S(λ) jm(Ω) are orthogonal on a unit sphere, i.e., /integraldisplay ΩS(λ) jm(Ω)·S(λ/prime) j/primem/prime(Ω)dΩ =Njλδjj/primeδmm/primeδλλ/prime, (A 8) with the norm Njλ= (1−δ−1λ) Πj+δ−1λ, (A 9) Πj=j(j+ 1). (A 10) Since the Maxwell equations do not admit the existence of magnetic monopoles, the zeroth order spherical harmonic functions Y00(Ω) and S(λ) 00(Ω) are excluded from the parameterizations of the scalar magnetic potential (2.10–2.12) and magnetic induction vector (2.19), respectively. From definitions (A 4–A 6) the following relations can be easily derived, S(λ) j−m(Ω) = ( −1)mS(λ) jm(Ω), (A 11) er×S(0) jm(Ω) = −S(1) jm(Ω), (A 12) er×S(−1) jm(Ω) = 0, (A 13) er×S(1) jm(Ω) = S(0) jm(Ω), (A 14) grad [f(r)Yjm(Ω)] =df(r) drS(−1) jm(Ω) +f(r) rS(1) jm(Ω), (A 15) curl/bracketleftBig f(r)S(0) jm(Ω)/bracketrightBig =−Πjf(r) rS(−1) jm(Ω)−/parenleftbiggd dr+1 r/parenrightbigg f(r)S(1) jm(Ω),(A 16) curl/bracketleftBig f(r)S(−1) jm(Ω)/bracketrightBig =−f(r) rS(0) jm(Ω), (A 17) curl/bracketleftBig f(r)S(1) jm(Ω)/bracketrightBig =/parenleftbiggd dr+1 r/parenrightbigg f(r)S(0) jm(Ω), (A 18) 75 div/bracketleftBig f(r)S(0) jm(Ω)/bracketrightBig = 0, (A 19) div/bracketleftBig f(r)S(−1) jm(Ω)/bracketrightBig =/parenleftbiggd dr+2 r/parenrightbigg f(r)Yjm(Ω), (A 20) div/bracketleftBig f(r)S(1) jm(Ω)/bracketrightBig =−Πjf(r) rYjm(Ω). (A 21) Byf(r) we mean an arbitrary continuous function of radius r. 76 A. SPHERICAL HARMONIC FUNCTIONS 77 Appendix B Existence and uniqueness of the integral solution Here we prove the existence and uniqueness of the time discretized, weak formulation of the electromagnetic induction problem (2.51–2.53). The boundary conditions imposed on the external part of the field at the surface are not common in the variational formu- lations of electromagnetic problems, since they cannot be expressed without the use of the spherical harmonic parameterization. Therefore first we reformulate the problem as a standard Dirichlet boundary value problem (DBVP) and then we show that it is equiv- alent to the formulation (2.51–2.53) which is more suitable for the realistic configuration of the Earth and its surroundings. Findi+1b∈Hcurl,0∀i= 0,1,..., such that µ0 ∆t/parenleftbigi+1b,δB/parenrightbig +a0(i+1b,δB) =µ0 ∆t/parenleftBig iB−i+1ˆB,δB/parenrightBig −a1(iB,δB)−a0(i+1ˆB,δB) ∀δB∈Hcurl,0, (B 1) wherei+1ˆBis an arbitrary function from Hcurl, such that er×i+1ˆB=−er×gradU(ti+1) on∂G. (B 2) Theni+1B=i+1ˆB+i+1bis the solution of the time discretized DBVP. Note that the solutioni+1Bis separated in two parts, the functioni+1ˆBis constructed from the horizontal component of the total field at the surface and a homogeneous Dirich- let boundary condition is imposed oni+1b. Unlike in the formulation (2.51–2.53), the unknown functioni+1band the test functions δBare from the same functional space Hcurl,0in DBVP. It is elementary to prove that the DBVP satisfies the prerequisites of the Lax-Milgram theorem (Kˇ r´ ıˇ zek & Neittaanm¨ aki, 1990), i.e. that the form on the left hand side of (B 1) is sesquilinear, continuous, and Hcurl-elliptic, and the right hand side is a linear functional with respect to δB. Reminding the definitions of Hcurl-scalar product ( ·,·)Hcurland norm 78 B. EXISTENCE AND UNIQUENESS OF THE INTEGRAL SOLUTION || · || Hcurlfrom F.2, we show the continuity, /vextendsingle/vextendsingle/vextendsingleµ0 ∆t/parenleftbigi+1b,δB/parenrightbig +a0(i+1b,δB)/vextendsingle/vextendsingle/vextendsingle≤α1/vextendsingle/vextendsingle/vextendsingle/parenleftbigi+1b,δB/parenrightbig Hcurl/vextendsingle/vextendsingle/vextendsingle≤α1||i+1b||Hcurl||δB||Hcurl,(B 3) and theHcurl-ellipticity of the left hand side. µ0 ∆t/parenleftbigi+1b,i+1b/parenrightbig +a0(i+1b,i+1b)≥α2||i+1b||2 Hcurl, (B 4) Positive constants are denoted by αk. Here we have to note that since ρ0is spherically symmetric, it is possible to prove the ellipticity of the form with respect to the norm of the Sobolev space W1,2(G)3. This is possible even in the case of azimuthally symmetric resistivity (Martinec, 1997), but not generally. Continuity of the functional on the right hand side of equation (B 1) is evident from the inequality /vextendsingle/vextendsingle/vextendsingleµ0 ∆t/parenleftBig iB−i+1ˆB,δB/parenrightBig −a1(iB,δB)− −a0(i+1ˆB,δB)/vextendsingle/vextendsingle/vextendsingle≤/parenleftBig α3||iB||Hcurl+α4||i+1ˆB||Hcurl/parenrightBig ||δB||Hcurl.(B 5) Since the choice ofi+1ˆBis not unique, it remains to prove that the solution of (B 1) does not depend oni+1ˆB. Leti+1ˆB1andi+1ˆB2be two different functions from Hcurl, satisfying the same boundary condition (B 2). Leti+1b1andi+1b2be the solutions of (B 1), corresponding toi+1ˆB1andi+1ˆB2, respectively, and Bl=i+1ˆBl+i+1blbe the solutions of DBVP for l= 1,2. Then we can use a particular test function, δB=i+1B2−i+1B1, δB∈Hcurl,0. (B 6) Substituting (B 6) into (B 1) and subtracting the two equations (B 1) for l= 2,1 yields, µ0 ∆t/parenleftbigi+1B2−i+1B1,i+1B2−i+1B1/parenrightbig +a0(i+1B2−i+1B1,i+1B2−i+1B1) = 0.(B 7) Because of the continuity of the form µ0/∆t(·,·) +a0(·,·) proved above, this can be satisfied if, and only if, ||i+1B2−i+1B1||Hcurl= 0 (B 8) which implies, i+1B1=i+1B2. (B 9) We have proved that foriBknown from the previous time-step and for the boundary condition er×gradU(ti+1) there exists a unique solutioni+1Bof the DBVP (B 1). The Dirichlet boundary condition (B 2) is represented in the spherical harmonic pa- rameterization by equations (2.21) and (2.23) which prescribe the horizontal component 79 ofBon∂G. The vertical component of the magnetic induction vector on ∂G, as de- scribed by equation (2.22), is obtained as a part of the unique solution of the DBVP. The boundary conditions (2.24–2.25) have been introduced as a linear combination of (2.22–2.23). Therefore, the Galerkin system of linear equations (2.103) corresponding to problem (2.51–2.53) and the Galerkin system corresponding to the DBVP (B 1–B 2), and constructed by means of the spherical harmonic-finite element parameterization, are equivalent. One can be derived from the other by replacing particular rows by their linear combinations. The formulation (2.51–2.53) which prescribes the potential of the external field at the boundary, is equivalent to the Dirichlet formulation (B 1–B 2). 80 B. EXISTENCE AND UNIQUENESS OF THE INTEGRAL SOLUTION 81 Appendix C Tetrahedral mesh The time-domain, 3-D finite-element EM induction solver that was introduced in Chapter 3 is to a certain extent independent on the choice of the discretization of the computational domain. Any conforming tetrahedral mesh which discretizes the Earth’s mantle Gand its insulating surroundings Acan be employed, providing that it also approximates the three spherical surfaces, i.e., the core-mantle boundary, the mantle-atmosphere boundary, and the outer boundary, with sufficient accuracy. However, the quality of the mesh has great influence on the stability and convergence of the iterative matrix inversion schemes and on the accuracy of the solution. Everett (1997) suggested a high-quality tetrahedral mesh generator designed in particular for geophysical problems in spherical geometry. Here we briefly summarize his method and introduce two additional features that we employ in the presented validation of the FE EM induction solver. Firstly, the core is removed from the computational domain and the core-mantle boundary is created and suitably approximated. Secondly, local refinement which increases the resolution of the method in areas of particular interest is introduced. Two different algorithms are used to generate the interior mesh that discretizes G, aCMB≤r≤aand the exterior mesh that discretizes A,a≤r≤aO. The interior mesh generator is based on the algorithm of recursive subdivision of a tetrahedron ABCD into eight subtetrahedra (Liu & Joe, 1996) which is shown in the left panel of Figure C 1. New vertices,E,F,G,H,I, andJare added in the midpoints of the edges AB,AC,AD,BC, BD, andCD, respectively. Four subtetrahedra, AEFG ,BEHI ,CFHJ , andDGIJ are created at the corners. The remaining central octahedron EFGHIJ is then divided into four subtetrahedra by addition of one edge. From the three possible choices, EJ,FI, andGHsuch edge is selected that yields the best quality of resulting subtetrahedra. The quality is measured by a coefficient Q (Liu & Joe, 1996), defined as Q=12(3V)2 3/summationtext i/negationslash=jL2 ij, (C 1) whereVis the volume of tetrahedron and Lijare the lengths of its edges. Note that 0<Q≤1 andQ= 1 for a regular tetrahedron. One additional modification is applied to this algorithm in the spherical geometry. 82 C. TETRAHEDRAL MESH EA BD FC HG IJ Figure C 1: Left: Subdivision of a tetrahedron into eight subtetrahedra used in the generation of the interior mesh. Right: Initial approximation of the Earth’s mantle and application of the first level of subdivision. The innermost tetrahedra shown by pink color are removed from the discretized domain. Instead of placing the new vertex, let us say E, exactly in the middle of edge AB, rE=1 2(rA+rB), (C 2) we translate the new node along radius so that its radial distance is an average of radial distances of its parents. In other words, we replace equation (C 2) by rE=1 2(|rA|+|rB|)(rA+rB) |rA+rB|. (C 3) This assures that all nodes of the resulting interior mesh are placed on concentric spherical surfaces, and the inner and outer boundaries of the discretized spherical shell are well approximated. The discretization of the Earth’s mantle is started from eight tetrahedra. There are seven vertices in this initial state, one is placed in the centre, one at each pole, and four vertices are regularly distributed along the equator. When the first level of recursive subdivision is applied, as shown in the right panel of Figure C 1, eight innermost subtetrahedra representing the core are removed and the inner boundary is created. Since the radius of the Earth’s core is only slightly greater than one half of the Earth’s radius, aCMB.= 0.546a, (C 4) the vertices on the inner boundary can be placed at the radial distance aCMB without impairing the mesh quality considerably. Then the subdivision algorithm is recursively 83 F DE A BC Figure C 2: Subdivision of a triangular prism into three subtetrahedra used in the generation of the exterior mesh. applied, until required spatial resolution is reached. For the level of recursion l, there are Ni nnodes andNi ttetrahedra of the interior mesh, Ni n=1 6/parenleftbig 7 8l+ 15 4l+ 8 2l/parenrightbig + 2, (C 5) Ni t= 7 8l, (C 6) of whichNc nnodes are on the interior surface ∂G1andNc ttetrahedra have one triangular facet on∂G1, Nc n= 4l+ 2, (C 7) Nc t= 2 4l. (C 8) The mantle-atmosphere boundary ∂G2≡∂A1is approximated by Ns nnodes and Ns t triangles, Ns n= 4l+1+ 2, (C 9) Ns t= 2 4l+1. (C 10) The exterior mesh consists of nconcentric spherical shells with interfaces placed at a=r0<r1<...<r n=aO. The discretization of each shell is based on the triangulation of its inner surface. The outer surface is created by projecting the vertices from the inner surface (A,B, andCin Figure C 2) along radius to the outer surface ( D,E, andF), r{D,E,F }=rk+1 rkr{A,B,C }. (C 11) 84 C. TETRAHEDRAL MESH 05101520f(%) 0.00.20.40.60.81.0 Q 05101520f(%) 0.00.20.40.60.81.0 Q Figure C 3: Quality of the interior (left panel) and exterior (right panel) mesh. Histogram of Q defined by (C 1) is shown for regular mesh M5,20(see Table C.1). Addition of edges in radial direction divides the shell into triangular prisms ( ABCDEF ). By creating an additional edge at each tetragonal facet of the prism, a subdivision into three tetrahedra is obtained. This algorithm is started from the triangulation of the innermost boundary at awhich is taken from the discretization of the interior mesh introduced above. Since the number of tetrahedra in each shell of the exterior mesh is the same, the thickness of the shell should increase with radius in order to keep good tetrahedra quality throughout the domain. We use a recursive formula rk−rk−1=frk−1π 2l+1. (C 12) It defines the thickness of the shell which is also the length of the edges in the radial direction, as an f-multiple of the length of the horizontal edge on the equator. Here lis the refinement level of the interior mesh. One can immediately see that the exterior mesh comprises of Ne nnodes andNe ttetra- hedra, Ne n= (n+ 1)Ns n, (C 13) Ne t= 3nNs t, (C 14) includingNs nnodes andNs ttetrahedra with a boundary facet on ∂A1, andNo nnodes and No ttetrahedra with a boundary facet on ∂A2, No n=Ns n, (C 15) No t=Ns t. (C 16) Figure C 3 shows the histograms of the quality parameter Qof the interior and exterior meshes for the highest resolution that is used in the examples presented in Chapter 3. 85 Figure C 4: Regular interior (red) and exterior (blue) meshes at various resolution. Upper left: l= 1, n= 2. Upper right: l= 2, n= 3. Bottom left: l= 3, n= 6. Bottom right: l= 4, n= 10. The main advantage of the finite-element methods over spectral methods is the pos- sibility of local refinement that increases the resolution in areas of particular interest. Our implementation is based on the same algorithm by Liu & Joe (1996) that is used in the generation of the interior mesh. Firstly, we select an area where increased resolution is required and for each node of the mesh we decide whether it is inside or outside the selected area. Tetrahedra with all four vertices located inside the refinement zone are divided into eight subtetrahedra (Figure C 1). In order to keep the mesh conforming, particular attention is required at the boundary of the refinement area. Tetrahedra with 86 C. TETRAHEDRAL MESH EA BD FC G EA BD C Figure C 5: Local refinement at the boundary of the refinement area. Left: Vertices A,B,Care inside, four new tetrahedra are created. Right: Only vertices AandBare inside, the original tetrahedron is split in two. only three and two vertices in the refinement area are divided into four and two subte- trahedra, respectively, as is shown in Figure C 5. Tetrahedra with one, or zero vertices inside are not refined. Again, all new nodes added during the refinement procedure are created using the formula (C 3), rather than at the edges’ midpoints (C 2). Although the local refinement can be applied repeatedly, the quality of the tetrahedra at the boundary of the refinement zone necessarily deteriorates, as can be seen from comparison of Q histograms in Figures C 3 and C 6. This represents a trade-off between a well conditioned problem with large dimension and a low dimension problem, with an ill-conditioned matrix that can take long time, or even fail, to invert. 87 Meshl n f a O Local refinement dx M5,205 20 1.5 4.068a NO 241192 M5,335 33 1.5 10.245a NO 294466 (3500 km ≤r≤6871 km) ∧ Mr 4,104 10 1.5 3.816a ∧(10◦≤ϑ≤70◦)∧ 36300 ∧(5◦≤ϕ≤65◦) (3500 km ≤r≤6871 km) ∧ Mr 5,205 20 1.5 4.068a ∧(10◦≤ϑ≤70◦)∧ 290398 ∧(5◦≤ϕ≤65◦) Table C.1: Overview of meshes used in validation examples in Chapter 3. By l,n,f,a0, and dxwe respectively denote the recursion level used in the generation of the interior mesh, number of layers of the exterior mesh, thickening factor of the exterior mesh, resulting radius of the outer boundary ∂A2(in multiples of the Earth’s radius), and the resulting dimension of the FE matrix A(3.61). The intervals in the “Local refinement” column define the neighbourhood of the spherical inclusion, where one additional level of local refinement was applied. The refinement zone spans over parts of both the interior and exterior meshes. 05101520f(%) 0.00.20.40.60.81.0 Q 05101520f(%) 0.00.20.40.60.81.0 Q Figure C 6: Quality of the interior (left panel) and exterior (right panel) mesh. Histogram of Q defined by (C 1) is shown for a locally refined mesh Mr 5,20(see Table C.1). 88 C. TETRAHEDRAL MESH Figure C 7: Interior part of the mesh Mr 5,20. Note the local refinement applied in the vicinity of the nested spherical inclusion. 89 Appendix D Nodal finite element functions In this appendix we define the 3-D piecewise linear finite elements defined on a tetrahedral grid that are used in Chapter 3 to construct a finite-dimensional approximation of Sobolev spacesW1,2,h 01(G)k,W1,2,h 02(A), andW1,2,h(A) (see Table F.2). We also introduce particular formulae for volume and surface integrals of finite element products. We start with the definition of nodal finite element functions. Let rq= (xq,yq,zq) be theq-th node of a tetrahedral mesh. We assign it a finite element Ψ qdefined in each tetrahedron as a linear function of x,y, andzthat satisfies Ψq(xp,yp,zp) =δpq. (D 1) One can immediately see that Ψ qis non-zero only in those tetrahedra that include the q-th vertex, and that in each tetrahedron there are only four non-zero finite elements, each assigned to one of its vertices. In order to obtain particular forms of finite element functions and their first deriva- tives, it is useful to introduce the reference tetrahedron (Kˇ r´ ıˇ zek & Neittaanm¨ aki, 1990). Letr/prime= (η,ξ,ζ ) be the position vector in a local Cartesian coordinate system. In this coordinate system we define the reference tetrahedron KLMN by positioning its vertices at r/prime k= (0,0,0), r/prime l= (1,0,0), r/prime m= (0,1,0), r/prime n= (0,0,1), (D 2) as is shown in Figure D 1. For each tetrahedron Tκ, given by vertices rk,rl,rm, and rn(k/negationslash=l/negationslash=m/negationslash=n), we can define an affine mapping from the local coordinates to the global ones, r=rk+Rκ·r/prime, (D 3) where the matrix Rκis given by Rκ= (rl−rk,rm−rk,rn−rk). (D 4) 90 D. NODAL FINITE ELEMENT FUNCTIONS 0 ξ1ζ ηLMN K 11 Figure D 1: Reference tetrahedron in local coordinates ( η, ξ, ζ ). To each vertex of the reference tetrahedron we assign a finite element function ac- cording to D 1, Ψ/prime k(r/prime) = 1 −η−ξ−ζ, Ψ/prime l(r/prime) =η, Ψ/prime m(r/prime) =ξ, Ψ/prime n(r/prime) =ζ. (D 5) The first order derivatives of the reference finite elements expressed in the local coordinate system are, grad/primeΨ/prime k(r/prime) = ( −1,−1,−1), grad/primeΨ/prime l(r/prime) = (1,0,0), grad/primeΨ/prime m(r/prime) = (0,1,0), grad/primeΨ/prime n(r/prime) = (0,0,1). (D 6) Since the transformation (D 3) is linear, we can write Ψq(r) = Ψ/prime q(r/prime), (D 7) grad Ψ q(r) =R−T κ·grad/primeΨ/prime q(r/prime), (D 8) whereq∈ {k,l,m,n }andR−T κdenotes the inverse transpose matrix to Rκ. Now we can easily express the products of the finite elements and their derivatives, integrated over the tetrahedron. Namely, for each tetrahedron Tκand for each pair of 91 finite elements Ψ p, Ψqthat are non-zero on Tκ, we define Ipq κ=/integraldisplay TκΨpΨqdV=|detRκ|/integraldisplay T/primeκΨ/prime pΨ/prime qdV/prime= =1 120|detRκ|(1 +δpq) for rp,rq∈Tκ, (D 9) Jpq κ=/integraldisplay TκΨpgrad Ψ qdV=|detRκ|R−T κ·grad/primeΨ/prime q/integraldisplay T/primeκΨ/prime pdV/prime= =1 24|detRκ|R−T κ·grad/primeΨ/prime qforrp,rq∈Tκ, (D 10) Kpq κ=/integraldisplay Tκgrad Ψ p·grad Ψ qdV= =|detRκ|(R−T κ·grad/primeΨ/prime p)·(R−T κ·grad/primeΨ/prime q)/integraldisplay T/primeκdV/prime= =1 6|detRκ|(R−T κ·grad/primeΨ/prime p)·(R−T κ·grad/primeΨ/prime q) for rp,rq∈Tκ.(D 11) We will also express the surface integrals of finite elements over the mantle-atmosphere interface∂G2≡∂A1. The interface is approximated by triangular facets Fη⊂Tκ. For each facetFηwe define the unit normal nηpointing approximately in the radial direction, nη≈eronFη(note that eris not constant on Fη). Thus we approximate the outer normals n≈nηon∂G2andn≈ −nηon∂A1. Then we can introduce integrals Lpq η=/integraldisplay FηΨpnη×grad Ψ qdS= =1 3|Fη|nη×/parenleftbig R−T κ·grad/primeΨ/prime q/parenrightbig forrp,rq∈Fη, (D 12) where |Fη|is the surface of facet Fη. Although grad Ψ qis non-zero on Fηalso for the finite element corresponding to node rq∈Tκ\Fη, it is parallel to nηsince Ψ q= 0 onFη. Note that ex·Lpq η=/integraldisplay FηΨpnη·curl (Ψ qex)dS, (D 13) and similarly for ey·Lpq ηandez·Lpq η. This allows to express the surface integrals in both equations (3.48) and (3.50) by means of Lpq η. Finally, it is possible to show that Lpq η1+Lqp η1=−/parenleftbig Lpq η2+Lqp η2/parenrightbig , (D 14) 92 D. NODAL FINITE ELEMENT FUNCTIONS for any neighbouring facets Fη1,Fη2with common vertices rp,rq∈Fη1∩Fη2,p/negationslash=q, since both expression depend only on rq−rp. Due to this relation the surface integrals do not violate symmetry of matrices A,B(3.63–3.64). 93 Appendix E Moving least squares interpolation The solution of the 3-D finite element formulation of the EM induction problem intro- duced in Chapter 3 is obtained in the form of electromagnetic potentials A, Φ, andU. To express the solution by means of a quantity measurable by ground based and satel- lite born magnetometers, the magnetic induction vector B, the first order derivatives of potentials, namely curl Aand gradUare needed. The potentials are parameterized by piecewise linear base functions with piecewise constant analytical derivatives. Therefore, the magnetic induction Bobtained by direct differentiation of the base functions would change discontinuously across the boundaries of neighbouring tetrahedra. We prefer in- stead the moving least squares method (MLS) which is used for numerical differentiation of irregularly sampled fields in fluid dynamics (Gossier, 2001). The MLS method is based on approximation of a scalar function f(x,y,z ) in the vicinity of a point ( x0,y0,z0) by a polynomial ˜fofj-th order, f(x,y,z )≈˜f(x,y,z ) =k+l+m=j/summationdisplay k,l,m ≥0fklm/parenleftbiggx−x0 a/parenrightbiggk/parenleftbiggy−y0 a/parenrightbiggl/parenleftbiggz−z0 a/parenrightbiggm , (E 1) where the Earth’s radius is a natural choice for the scaling factor a. Let gridpoints (xi,yi,zi),i= 1,...,n be thennearest neighbours of point ( x0,y0,z0). The values fi=f(xi,yi,zi) are known. Then we determine the coefficients fklmof the polynomial (E 1) by minimizing the least squares misfit χ2=n/summationdisplay i=1wi/parenleftBig fi−˜f(xi,yi,zi)/parenrightBig2 = min, (E 2) where the importance of each node can be optionally weighted by wi, usually a decreasing function of distance from ( x0,y0,z0). Problem (E 1–E 2) can be easily arranged in a matrix form. The condition number of the least squares matrix depends on the number and position of neighbours and often the matrix is ill conditioned. Following Press et al. (1992) we use the singular value decomposition (SVD) method to obtain the minimum norm solution of (E 1–E 2). During this process, the non-significant data can be detected by their small singular values and excluded from the fit. Once the coefficients fklmare found, the first-order derivatives of 94 E. MOVING LEAST SQUARES INTERPOLATION -30-20-10010[nT] 400050006000 r [km]-30-20-10010[nT] 400050006000 r [km]-30-20-10010[nT] 400050006000 r [km]020406080100120 400050006000 r [km]020406080100120 400050006000 r [km]020406080100120 400050006000 r [km]024 400050006000 r [km]024 400050006000 r [km]024 400050006000 r [km] Figure E 1: The high-resolution FE solution at t= 72 h for nested spheres conductivity model (see Figure 3.5 for details) obtained by direct differentiation of base functions (shown by crosses) and by MLS fit of quadratic polynomials into 10 nearest neighbours (dots). Note the large oscillation of the Bϕ component (right plate) yielded by the former method. ˜fat (x0,y0,z0) are simply, ∂˜f ∂x/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle (x0,y0,z0)=f100, (E 3) ∂˜f ∂y/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle (x0,y0,z0)=f010, (E 4) ∂˜f ∂z/vextendsingle/vextendsingle/vextendsingle/vextendsingle/vextendsingle (x0,y0,z0)=f001. (E 5) Applying the MLS method on Ax,Ay,Azin the Earth’s mantle and on Uin the insulating surroundings yields the magnetic induction vector B. Detailed discussion of importance of various parameters on the MLS fit can be found in Gossier (2001). We experienced best results by using non-weighted ( wi= 1) fit of a quadratic polynomial (j= 2) inton= 10 neighbours. An example comparing the induction vector Bobtained by direct differentiation of base functions and by MLS fit is shown in figure E 1. 95 Appendix F Notation F.1 Used symbols rposition vector (r,ϑ,ϕ )spherical coordinates: radius ( r≥0), colatitude (0 ≤ϑ≤π), and longitude (0 ≤ϕ≤2π) Ωabbreviated notation of angular coordinates, Ω = ( ϑ,ϕ) (x,y,z )Cartesian coordinates er,eϑ,eϕunit base vectors of spherical coordinate system ex,ey,ezunit base vectors of Cartesian coordinate system Gcomputational domain filled with conductor (sphere or spherical shell) ∂G boundary of G; ifGis a spherical shell, then ∂G=∂G1∪∂G2, where ∂G1,∂G2is the inner and the outer boundary, respectively Aspherical shell surrounding Gwhich is filled with insulator ∂A boundary of A,∂A=∂A1∪∂A2, where∂A1, and∂A2is the inner and the outer boundary, respectively nouter normal to a boundary; in spherical geometry n=eron∂G2∪∂A2, while n=−eron∂G1∪∂A1 aEarth’s radius, a= 6371 km aCMB radius of the Core-Mantle boundary, aCMB.= 0.546a aOradius of the outermost boundary ∂A2 Yjm(Ω) scalar spherical harmonic functions S(λ) jm(Ω) vector spherical harmonic functions ψk(r)piecewise linear 1-D nodal finite elements Ψk(r)piecewise linear 3-D nodal finite elements defined on a tetrahedral grid 96 F. NOTATION χk(r)piecewise linear functions based on ψk(r) to express curl ψk(r)S(λ) jm(Ω) t, ∆ttime, time-step (s) σelectrical conductivity (S/m) ρelectrical resistivity (Ω m) µ0magnetic permeability, µ0= 4π10−7H/m B magnetic induction vector (T) H magnetic field intensity vector (A/m) Eelectrical field intensity vector (V/m) jelectrical current density vector (A/m2) j(e)electrical current density of primary inducing field of external origin (A/m2) jSelectrical current surface density vector (A/m) Avector magnetic potential (A) Uscalar magnetic potential (A) U(e),U(i)scalar magnetic potential of external (primary, inducing) and internal (secondary, induced) field φscalar electrical potential (V) Φtime-integrated scalar electrical potential (V s) F.2 Overview of functional spaces C(G)kspace of scalar ( k= 1) and vector ( k= 3) functions continuous on the domainG Cn(G)kspace of scalar ( k= 1) and vector ( k= 3) functions whose classical derivatives up to the n-th order belong to C(G)k C1(/angbracketleft0,∞))kspace of scalar ( k= 1) and vector ( k= 3) functions continuously differ- entiable with respect to ton interval /angbracketleft0,∞) D0(G)space of smooth (infinitely differentiable) scalar functions with compact support inGwhich are zero on ∂G L2(G)kspace of scalar ( k= 1), vector ( k= 3), and tensor ( k= 9) functions F.2. OVERVIEW OF FUNCTIONAL SPACES 97 that are square-integrable in G; the scalar products (f,g) =/integraldisplay GfgdV, (f,g) =/integraldisplay Gf·gdV, of functions from L2(G)1andL2(G)3, respectively, induce the norms ||f||L2= (f,f)1 2, ||f||L2= (f,f)1 2 L+ ∞(G)space of positive ( f > 0) scalar functions bounded in G; the norm in L+ ∞(G) is defined as ||f||L∞= max G|f| W1,2(G)kSobolev’s space of scalar ( k= 1) and vector ( k= 3) functions from L2(G)kwhose partial derivatives belong to L2(G)3k; the scalar products (f,g)W1,2=/integraldisplay G[fg+ gradf·gradg]dV, (f,g)W1,2=/integraldisplay G/bracketleftbig f·g+ grad f: (grad g)T/bracketrightbig dV, of functions from W1,2(G)1andW1,2(G)3, respectively, induce the norms ||f||W1,2= (f,f)1 2 W1,2, ||f||W1,2= (f,f)1 2 W1,2 W1,2 01(G)kspace of scalar ( k= 1) and vector ( k= 3) functions from W1,2(G)kthat are zero on the inner boundary ∂G1 W1,2 02(A)space of scalar functions from W1,2(A) that are zero on the outer bound- ary∂A2 Hcurl space of vector functions from L2(G)3whose rotation also belongs to L2(G)3; the scalar product (f,g)Hcurl=/integraldisplay G[f·g+ curl f·curlg]dV, induces the norm |f||Hcurl= (f,f)1 2 Hcurl Hcurl,0space of vector functions from Hcurlwhose tangential component is zero on the boundary ∂G 98 F. NOTATION