Multiscale and Upscaling Methods forGeophysical Electromagnetic Forward ModelingbyLuz Ange´lica Caudillo MataB.Sc. Computer Science, University of Guanajuato, 2004M.Sc. Computer Science and Industrial Mathematics,Mathematics Research Center (CIMAT), 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDoctor of PhilosophyinTHE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES(Geophysics)The University of British Columbia(Vancouver)October 2017c© Luz Ange´lica Caudillo Mata, 2017AbstractAccurate and efficient simulation of electromagnetic responses in realistic geo-physical settings is crucial to the exploration, imaging, and characterization ofburied natural resources, such as mineral and hydrocarbon deposits. However, inpractice, these simulations are computationally expensive. The geophysical set-tings consider highly heterogeneous media and features at multiple spatial scalesthat require a very large mesh to be accurately represented. This results in asystem of equations to be solved that often exceeds the limits of average com-puters. Thus, the key is to reduce the problem size but retain the accuracy of theelectromagnetic responses.Upscaling and multiscale techniques have been successfully applied to theproblem of simulating fluid flow through heterogeneous porous media, where theyare able to drastically reduce the size of the resulting fine-mesh system by castingit into a coarse-mesh system that is much cheaper to solve, while achieving alevel of accuracy similar to that obtained with conventional discretization schemes.Recognizing the success that such techniques have had in fluid flow applications,this dissertation extends their use for application to electromagnetic modeling.In this dissertation, two new parallel simulation methods for the quasi-staticMaxwell’s equations in the frequency domain are proposed: an upscaling frame-work for the electrical conductivity, and a multiscale finite volume with oversam-pling method. Both methods are combined with an adaptive mesh refinementtechnique (OcTree) to boost their computational performance. The performanceof these methods is demonstrated by using field-inspired and synthetic examplesthat include a large electrical conductivity contrast.iiThis investigation shows that both proposed methods are feasible to tacklegeophysical electromagnetic problems, where being able to reduce the size of theproblem can be particularly advantageous when extended domains are consid-ered or when the mesh must capture the spatial distribution of the media hetero-geneity outside the region where the electromagnetic responses are measured.Furthermore, both methods are new contributions to the literature in the field ofcomputational methods in geophysical electromagnetics. Finally, both methodsincrease the current predictive and analytic capabilities by making the simulationof electromagnetic responses in larger and more complex geophysical settingsmore feasible than currently is possible.iiiLay SummaryThe focus of this investigation is to study how to reduce the computational costof simulating the behavior of electromagnetic fields in complicated geophysicalsettings. This type of simulations are crucial to the exploration of buried naturalresources (mineral, groundwater and hydrocarbon deposits) using electromag-netic methods. However, the size of the computation they involve often exceedsthe limits of average computers. This investigation proposes two innovative math-ematical alternatives (an upscaling and a multiscale with oversampling method)that achieve a great reduction of the problem’s size and the simulation cost with-out sacrificing much accuracy. It also opens the door to use such alternatives tocreate more powerful computational environments capable of simulating electro-magnetic fields in larger and more complex geophysical settings than currentlyis possible. Furthermore, this study is the first one in the literature to demon-strate the practical use of the proposed alternatives in the context of geophysicalelectromagnetic problems.ivPrefaceAll of the work presented henceforth was conducted at the Geophysical InversionFacility in the Department of Earth, Ocean and Atmospheric Sciences at the Uni-versity of British Columbia (UBC), Point Grey campus. This research resulted intwo peer-reviewed publications and in three expanded conference proceedings. Ipresented this research in eleven international conferences.Chapter 3 presents extended and revised versions of the publications:1. L. A. Caudillo-Mata, E. Haber, L. J. Heagy, and C. Schwarzbach, A frame-work for the upscaling of the electrical conductivity in the quasi-static Maxwell’sequations, J. Comput. Appl. Math., vol. 317, pp. 388 to 402, 20172. L. A. Caudillo-Mata, E. Haber, L. J. Heagy, and D. W. Oldenburg, Numericalupscaling of electrical conductivity: A problem specific approach to gener-ate coarse-scale models, in SEG Technical Program Expanded Abstracts2014, 2014, pp. 680 to 684.I was the lead investigator, responsible for all major areas of concept formation,background investigation, software development, simulation process, analysis ofresults, as well as manuscript composition. D. W. Oldenburg and L. J. Heagywere involved in the early stages of concept formation, the design of the simula-tions presented in Section 3.5.1 and their corresponding analysis of results, andcontributed to manuscript edits for the SEG expanded abstract. L. J. Heagy alsocontributed to manuscript revision and edits for the peer-reviewed publication. C.Schwarzbach contributed to software development and the design of the simula-tions presented in Section 3.5. E. Haber was the supervisory author on this studyand contributed in all related processes to the project at all its stages.vChapter 4 presents extended and revised versions of the publications:1. L. A. Caudillo-Mata, E. Haber, and C. Schwarzbach, An oversampling tech-nique for the multiscale finite volume method to simulate electromagneticresponses in the frequency domain, Computational Geosciences, 2017,doi:10.1007/s10596-017-9647-y2. L. A. Caudillo-Mata, E. Haber, and C. Schwarzbach, An oversampling tech-nique for the multiscale finite volume method to simulate frequency-domainelectromagnetic responses, in SEG Technical Program Expanded Abstracts2016, 2016, pp. 981 to 985.3. L. A. Caudillo-Mata, E. Haber, and C. Schwarzbach. A multiscale finitevolume method with oversampling for geophysical electromagnetic simula-tions. In EAGE - ECMOR XV - 15th Eur. Conf. Math. Oil Recover. -Expand. Abstr., 2016. doi:10.3997/2214-4609.201601892.I was the lead investigator, responsible for all major areas of concept forma-tion, background investigation, simulation process, analysis of results, as well asmanuscript composition. C. Schwarzbach and I developed the necessary soft-ware to run the simulations; together we conceived and analyzed the experiments.I performed all the experiments presented in Section 4.4. E. Haber was the super-visory author on this study and contributed in all related processes to the projectat all its stages.viTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiGlossary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . 41.3 Aim and Scope of the Study . . . . . . . . . . . . . . . . . . . . 91.4 Overview of the Study . . . . . . . . . . . . . . . . . . . . . . . . 102 Mathematical Background . . . . . . . . . . . . . . . . . . . . . . . 122.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122.2 The Electromagnetic Forward Problem . . . . . . . . . . . . . . . 12vii2.3 The Quasi-static Maxwell’s Equations in the Frequency Domain . 162.4 The Mimetic Finite Volume Discretization . . . . . . . . . . . . . . 172.5 Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.6 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . 292.6.1 Upscaling Methods . . . . . . . . . . . . . . . . . . . . . 302.6.2 Multiscale Finite Element and Finite Volume Methods . . . 342.6.3 Links between Multiscale FE/FV and Algebraic Multigrid . 363 An Upscaling Framework for the Electrical Conductivity . . . . . . 403.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.2 A Least-Squares Formulation of the Upscaling Problem . . . . . . 413.3 Numerical Results in 1D . . . . . . . . . . . . . . . . . . . . . . 483.3.1 Simulations Using a Single Frequency . . . . . . . . . . . 483.3.2 Simulations Using Multiple Frequencies . . . . . . . . . . 523.4 A Local Upscaling Framework for 3D Simulations . . . . . . . . . 563.5 Numerical Results in 3D . . . . . . . . . . . . . . . . . . . . . . 653.5.1 Simulations on a Single Coarse-mesh Cell . . . . . . . . 653.5.2 Simulations Using the Synthetic Lalor Model . . . . . . . 713.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 793.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 804 A Multiscale Finite Volume Method with Oversampling . . . . . . . 824.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 824.2 An Overview of the Multiscale Finite Volume Method . . . . . . . 834.3 The Oversampling Method . . . . . . . . . . . . . . . . . . . . . 884.4 Numerical Results in 3D . . . . . . . . . . . . . . . . . . . . . . 924.4.1 Simulations for the Synthetic Lalor Model . . . . . . . . . 924.4.2 Simulations for Random Heterogeneous Isotropic Media . 1044.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1165 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . 1185.1 Summary of Work and Research Highlights . . . . . . . . . . . . 1185.1.1 Upscaling: Summary and Highlights . . . . . . . . . . . . 1225.1.2 Multiscale with Oversampling: Summary and Highlights . 123viii5.2 Impact to the Field of Computational Methods in Geophysical Elec-tromagnetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1245.3 Remaining Challenges and Future Research . . . . . . . . . . . 125Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129ixList of TablesTable 3.1 Magnitude of the magnetic field (H-field) datum and relative er-rors resulting from forward modeling using the fine-mesh andfour coarse-mesh conductivity models. Note: %∗ denotes per-centage of primary field. . . . . . . . . . . . . . . . . . . . . 53Table 3.2 Relative errors for four coarse-mesh conductivity models andfor five frequencies. . . . . . . . . . . . . . . . . . . . . . . . 54Table 3.3 Analytical expressions for the set of linearly independent bound-ary conditions used to locally generate data in the cuboid cellΩH,extk (one per edge). Note that~x = (x1,x2,x3) ∈ R3. . . . . . 59Table 3.4 Relative errors in Euclidean norm for the magnitude of the sec-ondary magnetic fluxes induced by the ground in the surveyarea, ∆~B. The relative errors are given in per cent. . . . . . . 76Table 4.1 Relative errors in Euclidean norm for the total, real and imag-inary parts of the secondary magnetic fluxes induced by themineral deposit in the survey area, ∆~Bdeposit. Note that pc standsfor padding cells. . . . . . . . . . . . . . . . . . . . . . . . . 102Table 4.2 Relative errors (per cent) in Euclidean norm for the magnitudeof the x, y and z components of the secondary magnetic fluxesinduced by the mineral deposit in the survey area, ∆~Bdeposit. Notethat pc stands for padding cells. . . . . . . . . . . . . . . . . 103xTable 4.3 Average run time in seconds per simulation required to com-pute ∆~Bdeposit on a two hexa-core Intel Xeon X5660 CPUs at 2.8Hz with 64 GB shared RAM using MATLAB. Setup time: timerequired to compute the local interpolation matrices and to as-semble the reduced system of equations to be solved. Solvetime: time to solve the reduced system of equations. Total time:the sum of the setup and solve times. . . . . . . . . . . . . . . 104Table 4.4 Relative errors in Euclidean norm for the total, real and imag-inary parts of the secondary magnetic fluxes induced by therandom heterogeneous conductive medium in the survey area,∆~Brandom. Note that pc stands for padding cells. . . . . . . . . . 113Table 4.5 Relative errors (per cent) in Euclidean norm for the magnitudeof the x, y and z components of the secondary magnetic fluxesinduced by the mineral deposit in the survey area, ∆~Brandom. Notethat pc stands for padding cells. . . . . . . . . . . . . . . . . 114Table 4.6 Average run time in seconds per simulation required to com-pute ∆~Brandom on a two hexa-core Intel Xeon X5660 CPUs at 2.8Hz with 64 GB shared RAM using MATLAB. Setup time: timerequired to compute the local interpolation matrices and to as-semble the reduced system of equations to be solved. Solvetime: time to solve the reduced system of equations. Total time:the sum of the setup and solve times. . . . . . . . . . . . . . . 115xiList of FiguresFigure 1.1 Sketch of a geophysical large-loop EM experiment used inmining exploration and the induced secondary EM fields gen-erated by the mineralized conductive body. Figure courtesy ofBoliden Group (www.boliden.com). . . . . . . . . . . . . . . 2Figure 1.2 Geophysical inversion of EM data for natural resource explo-ration. The output is an electrical conductivity model of theunderground structures that were surveyed. Figure taken fromthe EM geosci website: https://em.geosci.xyz, which distributesits content under Creative Commons 4.0. . . . . . . . . . . . 3Figure 1.3 Example of a realistic geophysical EM setting that considers alarge computational domain, features that vary at multiple spa-tial scales (millimeters, meters, and kilometers), and a widevariation over several orders of magnitude of the physical prop-erties of the medium. Figure not drawn to scale. Figure cour-tesy of Rockware, Inc (www.rockware.com). . . . . . . . . . 6Figure 2.1 Control volume cell showing the staggered discretization for ~Eand ~Js on its edges, ~B on its faces, and the medium parametersµ and Σ at the cell centers. . . . . . . . . . . . . . . . . . . 20Figure 2.2 Schematic representation of the nearest-neighbor interpola-tion process for approximating (a) the edge inner product (Σ~E, ~W )and (b) the face inner product (µ−1~B,~F) in the cuboid cellΩhi, j,k. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23xiiFigure 3.1 Upscaling process. Maxwell’s equations are solved over thedomain Ω, where the electrical conductivity can vary over mul-tiple spatial scales and several orders of magnitude. The small-est spatial scale over which the conductivity varies defines thesize of the fine mesh on which the model is discretized. Thegoal of an upscaling procedure is to homogenize the conduc-tivity inside a subregion Ωup of Ω in order to construct an up-scaled electrical conductivity model suitable for simulation ona much coarser mesh. . . . . . . . . . . . . . . . . . . . . . 42Figure 3.2 Induction resistivity log: AA-05-01-096-11w4-0 from the Mc-Murray/Wabiskaw Oil Sands deposit well log database. (a)Discrete fine-scale conductivity model. Each diamond repre-sents a conductivity value on a uniform fine mesh of thickness25 cm. Each straight line represents a coarse layer of a uni-form mesh of thickness 10 m. The setup considers 320 finelayers and 8 coarse layers. (b) Resulting coarse-mesh con-ductivity models after applying four upscaling procedures: 1Dnumerical upscaling procedure (red solid line), and arithmetic(green dot line), geometric (magenta dot dash line), and har-monic (blue dash line) averages. . . . . . . . . . . . . . . . 51Figure 3.3 (a) Coarse-mesh electrical conductivity models obtained byusing the proposed upscaling method for different frequen-cies. The setup considers 320 fine layers and 8 coarse lay-ers. (b) Magnitude of magnetic field for each frequency, in %of primary field (%∗), resulting from forward modeling usingthe fine-mesh electrical conductivity model (black solid line),the different coarse-mesh conductivity models displayed in (a)(blue dash line), and the coarse-mesh models produced byusing arithmetic (red plus dot line), geometric (gray circle dotline), and harmonic averages (green square dot line). . . . . 55xiiiFigure 3.4 Local upscaling procedure for 3D settings. Left: fine-meshelectrical conductivity model and example of nested meshessetup. Right: extended domain (ΩH,extk ) for a given coarse-mesh cell (ΩHk ) and resulting anisotropic upscaled electricalconductivity (ΣHk ). . . . . . . . . . . . . . . . . . . . . . . . 57Figure 3.5 Non-zero component of each vectorial basis function ~Φl (de-fined in Table 3.3) plotted on a unitary cube. . . . . . . . . . . 60Figure 3.6 Cross sections of the fine scale conductivity model for (a) theisolated block and (b) the sheet. Conductive bodies are dis-played in red. The resistive background is displayed in blue.The coarse-mesh cell, ΩHk , which we aim to upscale is out-lined in white. . . . . . . . . . . . . . . . . . . . . . . . . . . 66Figure 3.7 Upscaled conductivity results found using the b and j criteriafor the conductive block model. Note that the upscaled con-ductivity in this case can be defined by a scalar, as the matrixrecovered was diagonal, with all diagonal elements being equal. 68Figure 3.8 Upscaled conductivity results found using the j and b criteriafor the conductive sheet. Since the upscaled conductivity is adiagonal 3×3 SPD matrix, it can be described by two positivescalars: σ1 = σ2 and σ3. σ1 = σ2 is shown in plot (a), and σ3is shown in plot (b). . . . . . . . . . . . . . . . . . . . . . . 70Figure 3.9 Subsurface part of the synthetic electrical conductivity Lalormodel and large-loop EM survey setup. The model is dis-cretized on a fine OcTree mesh with 546,295 cells. The con-ductivity varies over five orders of magnitude throughout thewhole model. . . . . . . . . . . . . . . . . . . . . . . . . . . 72Figure 3.10 Subsurface part of the synthetic electrical conductivity Lalormodel and large-loop EM survey setup. The model is dis-cretized on a coarse OcTree mesh with 60,656 cells. The con-ductivity varies over eight orders of magnitude throughout thewhole model. . . . . . . . . . . . . . . . . . . . . . . . . . . 74xivFigure 3.11 Magnitude of the secondary magnetic flux induced by the groundin the survey area, ∆, for the large-loop EM survey configura-tion. First and second rows show the results for 1 and 20 Hz,respectively. The left-hand panel shows the reference solutioncomputed using the MFV method on the fine OcTree meshwith 546,295 cells. The right-hand panel shows the upscaledsolution computed using 8 padding cells on the coarse OcTreemesh with 60,656 cells. . . . . . . . . . . . . . . . . . . . . . 78Figure 4.1 Schematic representation of the procedure to implement theMSFV method. . . . . . . . . . . . . . . . . . . . . . . . . . 83Figure 4.2 Schematic representation of the two main steps to implementthe oversampling method. . . . . . . . . . . . . . . . . . . . 89Figure 4.3 Subsurface part of the synthetic electrical conductivity Lalormodel and large-loop EM survey setup. The model is dis-cretized on a fine OcTree mesh with 546,295 cells. The con-ductivity varies over five orders of magnitude throughout thewhole model. . . . . . . . . . . . . . . . . . . . . . . . . . . 93Figure 4.4 Euclidean norm of total, real and imaginary parts of the refer-ence (fine-mesh) solution for the Lalor conductivity model perfrequency. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95Figure 4.5 Subsurface part of the synthetic electrical conductivity Lalormodel and large-loop EM survey setup. The model is dis-cretized on a coarse OcTree mesh with 60,656 cells. The con-ductivity varies over eight orders of magnitude throughout thewhole model. . . . . . . . . . . . . . . . . . . . . . . . . . . 96xvFigure 4.6 Real part of the z component of the secondary magnetic fluxinduced by the mineral deposit in the survey area, ∆Bzdeposit,for our large-loop EM survey at 100 Hz. (a) reference so-lution computed using the MFV method on the fine OcTreemesh with 546,295 cells. (b), (c) and (d): results using theMSFV+O method with 2, 4 and 8 padding cells on the coarseOcTree mesh, respectively. (e): results using the MSFV (with-out oversampling) method on the coarse OcTree mesh with60,656 cells. (f), (g) and (h): results using MFV with the con-ductivity model homogenized using arithmetic, geometric andharmonic averaging on the coarse OcTree mesh, respectively.All results are shown in picoteslas (pT) and plotted using thesame color scale. . . . . . . . . . . . . . . . . . . . . . . . 98Figure 4.7 Imaginary part of the z component of the secondary magneticflux induced by the mineral deposit in the survey area, ∆Bzdeposit,for our large-loop EM survey at 100 Hz. (a) reference solutioncomputed using the MFV method on the fine OcTree meshwith 546,295 cells. (b), (c) and (d): results using the MSFV+Omethod with 2, 4 and 8 padding cells on the coarse OcTreemesh, respectively. (e): results using the MSFV (without over-sampling) method on the coarse OcTree mesh with 60,656cells. (f), (g) and (h): results using MFV with the conductivitymodel homogenized using arithmetic, geometric and harmonicaveraging on the coarse OcTree mesh, respectively. All resultsare shown in picoteslas (pT) and plotted using the same colorscale. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99Figure 4.8 Subsurface part of the random heterogeneous isotropic elec-trical conductivity model and large-loop EM survey setup. Themodel is discretized on a fine OcTree mesh with 546,295 cells.The conductivity varies over eight orders of magnitude through-out the whole model. . . . . . . . . . . . . . . . . . . . . . . 106xviFigure 4.9 Euclidean norm of total, real and imaginary parts of the ref-erence (fine-mesh) solution for the random conductivity modelper frequency. . . . . . . . . . . . . . . . . . . . . . . . . . . 107Figure 4.10 Subsurface part of the random heterogeneous isotropic elec-trical conductivity model and large-loop EM survey setup. Themodel is discretized on a coarse OcTree mesh with 60,656cells. The conductivity varies over eight orders of magnitudethroughout the whole model. . . . . . . . . . . . . . . . . . . 108Figure 4.11 Real part of the z component of the secondary magnetic fluxinduced by the random heterogeneous conductive medium inthe survey area, ∆Bzrandom, for our large-loop EM survey at 100Hz. (a) reference solution computed using the MFV method onthe fine OcTree mesh with 546,295 cells. (b), (c) and (d): re-sults using the MSFV+O method with 1, 2 and 4 padding cellson the coarse OcTree mesh, respectively. (e): results usingthe MSFV (without oversampling) method on the coarse Oc-Tree mesh with 60,656 cells. (f), (g) and (h): results using MFVwith the conductivity model homogenized using arithmetic, ge-ometric and harmonic averaging on the coarse OcTree mesh,respectively. All results are shown in picoteslas (pT) and plot-ted using the same color scale. . . . . . . . . . . . . . . . . 111xviiFigure 4.12 Imaginary part of the z component of the secondary mag-netic flux induced by the random heterogeneous conductivemedium in the survey area, ∆Bzrandom, for our large-loop EM sur-vey at 100 Hz. (a) reference solution computed using the MFVmethod on the fine OcTree mesh with 546,295 cells. (b), (c)and (d): results using the MSFV+O method with 1, 2 and 4padding cells on the coarse OcTree mesh, respectively. (e):results using the MSFV (without oversampling) method on thecoarse OcTree mesh with 60,656 cells. (f), (g) and (h): resultsusing MFV with the conductivity model homogenized usingarithmetic, geometric and harmonic averaging on the coarseOcTree mesh, respectively. All results are shown in picoteslas(pT) and plotted using the same color scale. . . . . . . . . . 112xviiiGlossaryDOF Degrees of FreedomEM ElectromagneticFD Finite Difference discretization methodFE Finite Element discretization methodFV Finite Volume discretization methodMFV Mimetic Finite Volume discretization methodMSFV Multiscale Finite Volume MethodMSFV+O Multiscale Finite Volume Method with OversamplingPDE Partial Differential EquationSPD Symmetric Positive DefiniteUBC University of British ColumbiaxixAcknowledgmentsThis dissertation represents a long-awaited dream, which to be honest, there were sometimes I thought I will not be able to realize. In this space, I want to acknowledge to thosewho have provided me with words of wisdom to overcome my mental barriers in thosedifficult moments, as well as to those who have contributed to make my PhD experienceoutstanding. Two pages are not enough to express my deep gratitude to you all!At the top of the list is my advisor, Dr. Eldad Haber. Eldad, you have earned mydeepest gratitude, admiration, and respect not only as a scientist, but also as a humanbeing. Thank you for providing me with the right balance of freedom and guidance duringmy PhD; it allowed me to recognize and nurture my own ‘scientific voice’. Thank you fortaking me as a student, for providing me with continuous support (emotional, scientific,financial) when I have needed it the most, for pointing me to truly interesting and excitingprojects to work on, and for continuously challenging me and helping me to grow.I also thank the members of the GIF research group. I was truly blessed to be sur-rounded by such an extraordinary, diverse, and well-rounded group of people. I benefiteda lot from the can-do culture of our group. The standards that you set for yourselvesin the quality of your work helped me to set a high standard for me to try to live up to.Special thanks to Dr. Doug Oldenburg for his scientific vision and for his efforts (includingChristmas dinners, dragon boating, and hikes) on maintaining a cohesive group through-out time. Thanks as well to Christoph, Jenn, Klara, Kris, Gudni, Dominique, Sarah, Lars,Eran, Lindsey, and Wenke for their insightful contributions to improve my work and for fun.I also thank Drs. D. Oldenburg and M. Friedlander for being part of my supervisorycommittee and for their valuable contributions to improve my dissertation; to Drs. U.Ascher, C. Greif, A. Adem, M. Lamoureux, B. Homsey, and B. Wetton for their wisdomwhen I needed good advice; and to Dr. F. Herrmann for inviting me to work on geophysicalinverse problems when I started at UBC, which let me to pursue research in this field.Many thanks for the funding provided to pursue my studies by UBC through a Four-Year fellowship, EOAS through a H. W. Matthews scholarship, Mitacs through its Acceler-xxate program, and the members of the industrial consortium that support the GIF group.I also thank my other colleagues and co-workers in EOAS, IAM, and PIMS for creatingsuch an pleasant work environment and for their administrative and technical assistance.Special thanks go to H. Wason, C. Rada, D. Herrera, M. Joyce, R. Situma, D. Feighan, A.Van Slyck, D. Ruiz, and T. Nguyen.During my PhD, I did two internships at Lawrence Livermore (LLNL) and LawrenceBerkeley (LBL) National Laboratories. I thank my mentors there, Drs. P. Vassilevski and X.S. Li, respectively, for providing me with such amazing opportunities to work with them. Ialso thank the friends I made there for the scientific discussions we had (at LLNL: Delyan,Chack, Umberto, Max, Andrew, and Manuel. At LBL: Pieter and Mathias). Special thanksgo Drs. D. Awargal, F. Perez, and A. Belkacem at LBL for their generous professionaladvice.My life in Vancouver has been enriched by great friends, who have provided me with aspace to relax, to laugh, to connect, and to recharge. Thank you Alexandra, Yuri & Denis,Gabriela, Magdalena, Guillermina, Keiko, Claudia & Angel, Nidia & David, Male & Javier,and Diana. I also thank to my other close friends around the world who have continue tomake the effort to keep our long-distance friendship going.Without the continuous encouragement from my CIMAT mentors, Drs. Adolfo Sa´nchezValenzuela, Ignacio Barradas, and M. Tapia; as well as from my good friends and Englishteachers, L. Maxwell, G. & L. Henderson, and L. Nash, I would have not pursue a PhDabroad. Thank you Adolfo and Ignacio for being the first ones who believed in me andmade me believe that I had the potential to pursue a career in Industrial Mathematics. Ialso thank Drs. J. Nocedal, L. Jodar, J. A. de la Pen˜a, and C. Castillo Chavez, whoseadvice was crucial to identify a great doctoral program abroad aligned to my interests.There are simply no words to express my utmost gratitude to my parents, RamonaMata Bonilla and Isauro Jesu´s Caudillo Herrera, for all they have done for granting meaccess to an education, and then to ensure I thrive while I pursue it. Thank you to both ofyou, my two brothers (Missael and Hiram), my two nieces (Michelle and Ximena), and mytwo nephews (Axel and Alan), for being a continuous source of strength and inspirationand for being present in the good and bad moments despite the distance. I also thank thesupport of my parents-in-law: Elı´as Antonio Huchim and Beatriz Novelo.Finally, immense thanks to my beloved husband: Elı´as de Jesu´s Huchim Novelo, forenabling me every day to achieve my dreams. I have only received love, support andencouragement from you since I met you. Thanks as well for your valuable contributionsto improve my work, and for the overall fulfilling experience of sharing our life together! Toyou, my dearest friend, my eternal gratitude and admiration. You have become my home.xxiDedicationA mi querida familia: Elias, Ramona, Isauro, Missael, Hiram,Michelle, Ximena, Alan, Axel.[To my beloved family: Elias, Ramona, Isauro, Missael, Hiram,Michelle, Ximena, Alan, Axel.]A mi entrenador de ajedrez: Marco Antonio Rodrı´guez Ortega.[To my chess coach: Marco Antonio Rodriguez Ortega.]xxiiChapter 1IntroductionEverything can be taken from a person but one thing: the last of thehuman freedoms – to choose one’s attitude in any given set ofcircumstances, to choose one’s own way. — Viktor Frankl1.1 MotivationSince their inception in the early 1900s, geophysical electromagnetic (EM) meth-ods have been used in industrial, geoscientific and engineering applications world-wide to detect, locate and characterize buried natural resources of economic sig-nificance. Currently, the two dominant commercial uses of geophysical EM meth-ods are in exploration and monitoring of mineral deposits and hydrocarbon (oil andgas) reserves (see [91, 139, 153, 166, 179] and references within). Other com-mercial uses include exploration of geothermal and groundwater deposits (e.g.[42, 155]), detection of buried unexploded ordnance (e.g. [154]), and monitoringenvironmental impact, such as carbon dioxide sequestration (e.g. [42]).Geophysical EM methods infer underground structures from measurementsof electric and/or magnetic fields that are acquired through a variety of surveysettings [144, 161]. A geophysical EM survey consists of one or several EM ex-periments. Figure 1.1 shows an example of a typical EM experiment setup usedin the mining exploration.A typical EM experiment starts by energizing the ground using an electro-magnetic source, which can be harmonic or transient, natural or controlled (cf.1[144]). The energy propagates throughout the subsurface and induces secondaryelectromagnetic fields in the regions where there is a (high) contrast in electricalconductivity (Figure 1.1). Electrical conductivity is the physical property of theground that quantifies the ability of a material to allow the flow of electric currents,and it is used in exploration geophysics to characterize Earth’s (conductive or re-sistive) materials [161]. Next, EM measurements or data containing informationof the distribution of the subsurface conductivity are recorded by placing receiversat the Earth’s surface, in the air, in the sea, or in wells. The experiment is typicallyrepeated for various source and receiver locations, resulting in a large volume ofdata. Subsequent processing, modeling and inversion of EM data is performed toconstruct detailed subsurface electrical conductivity models. The geological inter-pretation and assessment of these conductivity models allows the understandingof the spatial distribution of the underground structures, and it ultimately leadsto better-informed exploration decisions. Figure 1.2 shows the concept of a geo-physical inversion of EM data.Figure 1.1: Sketch of a geophysical large-loop EM experiment used in min-ing exploration and the induced secondary EM fields generated bythe mineralized conductive body. Figure courtesy of Boliden Group(www.boliden.com).2Figure 1.2: Geophysical inversion of EM data for natural resource explo-ration. The output is an electrical conductivity model of the under-ground structures that were surveyed. Figure taken from the EM geosciwebsite: https://em.geosci.xyz, which distributes its content under Cre-ative Commons 4.0.Modeling and inversion of geophysical EM data are powerful interpretationtools; however, using them as part of exploration programs in industry is both timeconsuming and computationally challenging [19, 68, 134, 139, 153, 166, 179].Inversion of EM data is achieved by solving a non-linear and (often) large-scale inverse problem using optimization techniques [138, 159]. The major bot-3tleneck of performing an inversion is to solve the forward problem. The EM for-ward problem is the mathematical procedure by which, given information of thephysical properties of the medium and the sources, we can predict (simulate) EMdata through numerically solving the corresponding governing Maxwell’s equa-tions [156, 159, 169]. That is, the solution of the forward problem involves anumerical process for solving a partial differential equation (PDE) or a system ofPDEs (depending on the formulation of Maxwell’s equations that is chosen to workwith). In the geosciences, solving a forward problem is also referred as forwardmodeling or just modeling. Besides playing a key role within an inversion, the EMforward problem is by itself crucial to simulate physical processes governed byEM induction phenomena — which has a rich spectrum of applications in differ-ent fields of knowledge, including, medical imaging and the various engineeringbranches (electrical, mechanical, materials, to name a few) [58, 78, 129, 156].Even though there has been significant progress in advancing existing com-putational EM capability, solving the forward problem for realistic, large-scale geo-physical settings remains a computational challenge (see [19, 68, 134, 139, 153,166, 179] and references within). As we will see in the next section, such chal-lenge consists on dealing effectively with the elevated computational cost (bothin CPU memory and run time) of this type of simulations. Therefore, working to-wards faster and more accurate computational solutions for the geophysical EMforward problem will have a tangible impact in the overall processes of imagingand monitoring buried natural resources using geophysical EM methods.1.2 Problem StatementAccurate and efficient simulation of EM responses — EM fields and fluxes —through solving the EM forward problem for realistic, large-scale geophysical set-tings can be a very computationally expensive task, specially when conductedwithin an EM inversion procedure (see [22, 68, 134, 179] and references within).The major challenge in practice to perform this type of (forward) simulation is thesize of the computation it involves.Realistic geophysical settings typically consider large computational domainsin 3D; geologic, topographic, and other simulation-related features that vary at4multiple spatial scales; and a wide variation over several orders of magnitude ofthe physical properties (e.g. electrical conductivity) of the heterogeneous media.Figure 1.3 shows one example of a complex geophysical setting that can be en-countered in practice, where some features in the model can vary from millimetersupward (e.g. the thickness of the casing of wells), while the simulation domain canbe on the order of hundreds of kilometers. In addition, the electrical conductivityof Earth’s materials varies over many orders of magnitude [161]. For the exampleshown in Figure 1.3, the values of the electrical conductivity can range from 10−8to 106 S/m, when the setting considers air and steel well casing, which correspondto the features in the model that are less and most conductive, respectively.Since all of the features considered in a geophysical setting can have a sig-nificant impact on the behavior of the EM responses of interest (that are typicallymeasured at the surface), if we need to obtain an accurate approximation to theresponses, the mesh used in traditional discretization techniques, such as finitevolume or finite element, should be able to capture the structure of the hetero-geneity present in the setting with sufficient detail. This need leads to the useof a very large mesh to discretize the model, which results into solving a verylarge, and often very ill-conditioned, system of equations — in some cases, inthe order of millions (106), or even larger than billions (109) of unknowns. For thegeophysical setting shown in Figure 1.3, the system of equations (depending onthe level of detail considered in the model) can lead to a system in the order of1012 of unknowns. Such a large system of equations require specialized comput-ing resources (e.g. clusters) to be solved. When an EM simulation is conductedin practice (e.g. for different frequencies or survey configurations) or within anEM inversion procedure (where there is the need to deal with adjoint operatorsper frequency and/or source), several forward EM simulations must be conducted(the book in [68] gives a more detail estimate of the number of linear systems tobe solved). This can lead to a very computational expensive process overall, andtherefore it is of interest to reduce effectively the computational cost of individualEM simulations.Conventional approaches used to reduce the computational cost of EM for-ward simulations are: parallel computing, adaptive mesh refinement, and iterativesolvers. The combination of these three approaches is one of the most efficient5Figure 1.3: Example of a realistic geophysical EM setting that considersa large computational domain, features that vary at multiple spatialscales (millimeters, meters, and kilometers), and a wide variationover several orders of magnitude of the physical properties of themedium. Figure not drawn to scale. Figure courtesy of Rockware,Inc (www.rockware.com).ways currently used to tackle this problem [68, 134]. However, the ever increasingneed to accurately simulate larger and more complex settings, which requires usto be able to include more detail in the modeling stage, limits their sustainability[22]. As we see next, each individual approach faces some issues to reduce thesize of the computation, in terms of both CPU memory and time requirements.Parallel computing drastically reduces the run time of the simulation by break-ing the problem into discrete pieces of work that can be carried out simultaneouslyusing multiple processors; however, the size of the system of equations is not re-duced. The workload is merely distributed among processors. See [100, 134] fora recent review of high-performance computational strategies for EM modeling.Adaptive mesh refinement (AMR) approaches have been used to overcome(reduce) the computational cost of EM modeling in geophysical applications (e.g.6[72, 86, 108, 118, 149]). They do so by refining selected areas of the mesh to meetcertain error bounds based on the requirements of the problem. AMR approacheshave produced accurate approximations to the EM responses at an affordablecost; however, when used to simulate geophysical EM responses in heteroge-neous settings, the mesh must still capture the relevant features in the spatialdistribution of the media heterogeneity both inside and outside the region wherewe measure the EM responses. This restricts the ability of AMR approaches toreduce the size of the system to be solved.Efficient iterative solvers and preconditioners have been developed for thelarge, sparse, and ill-conditioned system of equations that results from the dis-cretization of the Maxwell’s equations (e.g. [63, 64, 68, 70, 71, 83, 134, 172]).These solvers lead to substantial savings on time and memory usage, as they arebased on optimized (sparse) matrix-vector product multiplications and do not re-quire to store the matrix of the system [146]. Nevertheless, using iterative solversto solve the EM forward problem as part of an EM inversion routine suffer a ma-jor drawback: the system needs to be solved multiple times for the number ofsources and the different frequencies or time steps in the geophysical EM surveyconfiguration [68], which can be in the order of thousands or even millions (e.g.an airborne survey) [144]. In such cases, being able to use direct solvers to de-compose the matrix of the system of equations is the most cost-effective way toovercome this issue [43, 68]. However, a matrix factorization requires us to beable to store the matrix.The issues described suggest the need for a new computational frameworkcapable of reducing the size of the system of equations more effectively. Fur-thermore, the new framework needs to be able to compute accurate approxima-tions for the EM problem at a coarser resolution, such that the impact caused bysmall-scale features in the model is still included but without resolving for suchsmall-scale features. For instance, in the example provided in Figure 1.3, we cannot afford to mesh the model at the millimeter scale when the domain varies overhundreds of kilometers. Instead, we need an accurate approximation to the solu-tion at the meter scale that includes the effect due to the presence of the relevantmillimetric conductive feature (casing) in the model. Moreover, the new frame-work must be able to leverage most, if not all, of the advantages of the methods7described above. That is, it needs to run in parallel, be able to work with adaptivemesh refinement, and be accurate and cost-effective.The literature reports that upscaling and multiscale (finite volume and finiteelement) methods provide parallel methods to solve linear, elliptic PDE problemswith multiple scale features and highly heterogeneous settings. For such typeof PDE problems, both of these methods are able to compute large-scale physi-cal solutions that capture the small-scale effect accurately and efficiently withoutresolving for the small-scale features in the model [49, 51, 56, 88, 89].Favorably, upscaling and multiscale methods have been rigorously studied inthe field of petroleum engineering for the problem of modeling (single-phase) fluidflow through highly heterogeneous porous media, whose governing PDE modelis the steady-state, linear Poisson’s equation (cf. [6, 51]). Such a problem sharesseveral key challenges similar to the problem of simulating geophysical EM re-sponses in highly heterogeneous settings, namely, the governing PDE model inboth problems is linear, the simulation considers large-scale computational do-mains, features varying at multiple spatial scales, and a wide variation over severalorders of magnitude of the physical properties (e.g. permeability) of the media.For single-phase flow in porous media problems, these two methods have beensuccessfully used to drastically reduce the size of the fine-mesh system from thediscretization of the Poisson’s equation by constructing a coarse-mesh version ofthe system that is much cheaper to solve, while achieving a level of accuracy sim-ilar to that obtained with traditional discretization schemes (e.g. finite element orfinite volume) on a fine mesh [49–51, 56, 60, 80, 81, 88, 101, 122, 123, 132, 141,175].Since the use of upscaling and multiscale techniques has not been rigorouslyinvestigated for application to geophysical EM modeling, in this dissertation, I in-vestigate how these two techniques can be extended to do so. At present, de-veloping efficient multiscale methods using different discretization schemes andtailoring them for use to diverse applications is an active research area (e.g.[38, 39, 59, 87, 114] and references within).81.3 Aim and Scope of the StudyThe aim of this study is to investigate the applicability and feasibility of using up-scaling and multiscale (finite volume and finite element) techniques to efficientlysolve the geophysical EM forward problem for complex settings that include fea-tures varying at multiple spatial scales and several orders of magnitude. I developthe aim of this study by focusing on answering the following research questions:• How to extend the core mathematical ideas of some successful upscalingand multiscale methods developed for the problem of simulating fluid flow inhighly porous media for application to the problem of simulating geophysicalEM responses in highly heterogeneous conductive media? That is, how toextend some of the core mathematical procedures developed for a linearscalar PDE model (i.e. the Poisson’s equation) to a linear vector PDE model(i.e. the Maxwell’s equations)?• How do upscaling and multiscale methods perform for geophysical EM prob-lems with highly heterogeneous settings when they are combined with anadaptive mesh refinement technique? In particular, can these methods beused to drastically reduce the size of the system of equations to be solved?In order to answer these questions, the scope of this study is focused solelyon geophysical EM problems in the frequency domain where the quasi-static ap-proximation applies. The quasi-static approximation of EM fields refers to the factthat the contribution of the electric current displacement is negligible compared tothe electric current density when working at frequencies lower than 105 Hz [169].This is the case for a wide range of geophysical EM surveys used in practice.For example, we can model the EM responses of controlled-source EM surveys(e.g. large-loop surface and airborne surveys) and natural-source EM surveys(e.g. magnetotellurics and ZTEM surveys). The books in [144, 161] include athorough description for each of these surveys.Since most geophysical EM surveys in the quasi-static regime aim to charac-terize the electrical conductivity, we assume that the electrical conductivity is theprimary medium’s physical property of interest in this study [138, 144, 161, 179].91.4 Overview of the StudyGiven the background information and proposed problem statement in this intro-ductory chapter, the rest of this dissertation is organized as follows.Chapter 2 introduces the quasi-static Maxwell’s equations in the frequencydomain, which is the mathematical model I focus on in this dissertation; providesan overview of the mimetic finite volume discretization method, which is used asa building block to develop the upscaling and multiscale techniques proposed inthis work; and provides a summary of the relevant problem characteristics, theliterature review as well as related work.Chapter 3 develops an upscaling framework for the electrical conductivity ofthe quasi-static Maxwell’s equations in the frequency domain. The goal of suchframework is to construct accurate coarse-mesh conductivity models from givenfine-mesh ones that can be used to accurately simulate EM responses on acoarse mesh. The main contribution of this chapter is to pose upscaling as aparameter estimation problem to be solved on each coarse-mesh cell, which isfundamentally different than other upscaling formulations proposed in the litera-ture. This chapter starts by introducing the components for the framework thatare used to propose a general least-squares formulation for the upscaling prob-lem in geophysical EM applications. Then, I use a 1D example to illustrate thegeneral principle behind the framework and to show its performance to upscalewell log electrical conductivity data from the Canadian McMurray formation. Thisexample provides the perfect scenario to test our upscaling framework and tooffer a potential alternative to upscale well log electrical conductivity data in prac-tice. Afterwards, I adapt the upscaling formulation to construct coarse-mesh 3Danisotropic electrical conductivity models from given fine-mesh isotropic models.Finally, I show the performance of the adapted upscaling formulation for 3D prob-lems using two upscaling examples on a single coarse cell, and one example ofa synthetic electrical conductivity model from inversion results of field measure-ments over the Canadian Lalor mine. The last example also shows how the up-scaling method can be combined with OcTree (a type of adaptive mesh refinementtechnique) in a parallel environment to boost its performance.Chapter 4 develops a multiscale finite volume method with oversampling for10the quasi-static Maxwell’s equations in the frequency domain. The goal of thismethod is to provide more accurate approximations to the EM responses thanthose obtained with the multiscale finite volume technique (without oversampling)developed by [73] for tensor meshes and at a fraction of the cost of traditional dis-cretization techniques (such as mimetic finite volume) on a fine mesh. The maincontribution of this chapter is to show how the core mathematical ideas used todevelop one of the most successful oversampling techniques for flow in porousmedia applications can be extended to geophysical EM problems, as well as toshow how the multiscale finite volume method with and without oversampling canbe combined with OcTree meshes in a parallel environment to tackle more chal-lenging geophysical EM simulations. This chapter starts by providing a summaryof a multiscale finite volume method for geophysical EM problems. I discuss howthe accuracy of such method can be improved by complementing it with an over-sampling technique. Finally, I show the performance of the method using two 3Dsynthetic models: the Lalor model (introduced in Chapter 3), and one with randomisotropic heterogeneous media.Finally, Chapter 5 summarizes the work, discusses the contribution of thisinvestigation to the field of computational EM, and discusses the remaining chal-lenges and future research directions of the work presented in this dissertation.11Chapter 2Mathematical BackgroundLive as if you were to die tomorrow.Learn as if you were to live forever. — Mahatma Gandhi2.1 OverviewThis chapter provides a literature review of the numerical methods used to solvethe quasi-static EM forward problem in the frequency domain, introduces themathematical model I focus on in this dissertation, provides an overview of themimetic finite volume discretization method, discusses the existing solvers for theresulting system of equations obtained from the discretization of the Maxwell’sequations using the mimetic finite volume method, and concludes with a literaturereview of the related upscaling and multiscale techniques developed for flow inporous media problems that are used as a building block to develop the upscalingand multiscale with oversampling methods proposed in this doctoral dissertationfor application to geophysical EM problems.2.2 The Electromagnetic Forward ProblemThis section provides a literature review of the conventional methods that are usedto solve the geophysical quasi-static EM forward problem in the frequency domain.In Geophysics, the EM forward problem or forward modeling refers to themathematical procedure by which given information of the physical properties of12the medium (e.g. spatial distribution of both the electrical conductivity and mag-netic permeability) and the sources, we can predict (compute) geophysical EMdata (responses) by numerically solving the Maxwell’s equations [158, 169]. Theempirical Maxwell’s equations are widely accepted in the scientific community asthe mathematical model that governs all macroscopic EM phenomena [58, 156].Over the last four decades, the scientific community has been proposing effi-cient numerical methods to solve the EM forward problem for different geophys-ical applications [5, 7, 19, 34, 53, 68, 70, 76, 106, 142, 153, 157, 166, 179].The progress has happened on several fronts in the fields of scientific and high-performance computing. Nowadays, we have a broader set of numerical methods(e.g. adaptive mesh refinement and discretization techniques, solvers and pre-conditioners) to use as well as more powerful machines to run the simulations.This progress has enabled modelers to increase the complexity of the simula-tion. For example, we went from 1D settings to 3D settings, from small-scale tolarge-scale domains, from simple homogeneous or layered models that considertargets in the form of blocks to highly heterogeneous models that consider morecomplicated geological and topographical structures. In particular, performing 3Dsimulations has been possible due to the increase in computing power we havehad over the last two decades [134, 179]. Currently, we are still not at the pointwhere we can fully tackle computationally this problem. For example, runningan EM forward simulation on a geophysical setting as the one showed in Figure1.3 can not be handled without specialized computing resources (e.g. a clustersupercomputer), to which few people have access.To propose a better strategy to solve the geophysical EM forward problem, wefirst need to understand how it has been traditionally solved. The conventionalapproach to solve this problem consists of the following three major steps:1. Select a well-posed boundary value problem (in the sense of Hadamard)that represents the physical phenomenon of interest. A well-posed bound-ary value problem in the sense of Hadamard consist of a PDE model withboundary conditions that has a unique solution, which depends continu-ously on the input [77].The quasi-static Maxwell’s equations are the established mathematical PDE13model to simulate physical responses resulting from a geophysical EM ex-periment [68, 169]. This model can be formulated in the frequency or timedomain. For the sake of simplicity of exposition, I focus on frequency-domain EM problems. In addition, the model can be formulated in first-orderor second-order form. The literature reports that using the first-order formleads to more accurate results when both EM fields and fluxes need to becomputed [68]. In this dissertation, I work with the first-order formulation forthis reason.Section 2.3 provides full details about the boundary value problem I workwith in this dissertation.2. Discretize the quasi-static Maxwell’s equations using an appropriate meshand discretization method.The mesh used in this step can be structured, semi-structured or unstruc-tured. The literature reports that structured and semi-structured meshes areused more often than unstructured meshes for the geophysical EM forwardproblem. However, there is increasing research on how to use effectivelyunstructured meshes as well (e.g. [100, 115, 150] and references within).Structured and semi-structured meshes have the advantage of leading tomore accurate solutions, structured, symmetric and sparse matrices, aswell as being simpler to program and parallelize. Their main disadvantageis that they cannot model effectively complicated geometries (e.g. Figure1.3). The most common structured and semi-structured meshes used inpractice are: tensor meshes (i.e., any mesh that has a constant width alongthe entire axis such that it can be defined by a single width vector), or-thogonal locally refined meshes (e.g. QuadTree and OcTree), and logicallyorthogonal meshes (i.e., a distorted orthogonal mesh) [68].Conventional discretization methods for this step are: the finite differencediscretization method (FD) (e.g. [41, 99, 106, 178] and references within),the finite element discretization method (FE) (e.g. [100, 104, 115, 130, 131,149]), the finite volume discretization method (FV) (e.g. [54, 68, 100]), andthe Integral Equation method (e.g. [8, 126, 168]). The Integral Equation14method has been reported to be less competitive than the other methods[149]; it cannot deal with large conductivity contrasts neither with compli-cated geometries and it involves large dense matrices, which increase thecost of solving the system of equations. Currently, the most common dis-cretization methods used in practice for our problem are FD and FV meth-ods on orthogonal structured or semi-structured meshes, and FE on struc-tured and unstructured meshes.The state-of-the-Art is to use a mimetic discretization method as it resultson discrete solutions that preserve (mimic) the relevant underlying math-ematical and physical properties of the continuum PDE model on generalpolygonal and polyhedral meshes, and it leads to sparse and symmetricsystems of equations [68, 121]. Mimetic discretizations have been derivedfor FD, FE and FV methods [93–96, 104, 121, 131]. In this dissertation,we use the mimetic finite volume discretization method (covered in detailin Section 2.4) on tensor and OcTree meshes. Section 2.4 elaborates fur-ther on the advantages of using this discretization method as well as on thephysical properties that this mimetic discretization preserves for the quasi-static Maxwell’s equations.3. Solve the linear system of algebraic equations resulting from the previousstep to obtain an approximation to the EM responses. Section 2.5 discussesthe solver options for the system of equations that we actually solve in thisdissertation.Using the conventional approach outlined above to solve the quasi-static EMforward problem in practice is computationally very expensive. As discussed inSection 1.2, the cost comes from solving a very large system of equations thatresults from the discretization of a complex geophysical setting that includes fea-tures at multiple spatial scales (e.g. Figure 1.3), which often requires a very de-tailed mesh. As pointed out in [22, 23, 122], just building more powerful ma-chines does not constitute a sustainable methodology to solve this problem asthe amount of processing increases too steeply with the rise in problem size. Asdiscussed in Section 1.2, this study proposes two new methods, an upscalingmethod and a multiscale method with oversampling, to overcome this situation.152.3 The Quasi-static Maxwell’s Equations in theFrequency DomainThis section introduces the boundary value problem I focus on in this dissertation.More detailed discussions on EM theory for geophysical applications can be foundin [144, 156, 161, 169].Due to attenuation of the EM responses in the earth, geophysical EM sur-veys usually utilize frequencies in the range (0,105) Hz, in which case electricaldisplacement currents can be neglected [169]. For this scenario, the governingmathematical model is given by the first-order form of the quasi-static Maxwell’sequations ([68, 169]):∇× ~E + ıω~B = ~0, in Ω, (2.1)∇× µ−1~B−Σ~E = ~Js, in Ω, (2.2)where ~B denotes the magnetic flux density, ~E denotes the electric field intensity,~Js denotes the source term, ω denotes the angular frequency, ı is the unit imag-inary number, and Ω ⊂ R3 denotes the domain. The PDE coefficients, µ and Σ,denote the magnetic permeability and electrical conductivity, respectively. Notethat the system of equations (2.1)-(2.2) constitutes a linear, complex and vectorPDE system.The PDE coefficients, µ and Σ, can be scalars or Symmetric Positive Definite(SPD) tensors. For the sake of generality in the exposition, let us assume that bothcoefficients are 3×3 SPD matrices of the formµ(~x) =µ1(~x) µ4(~x) µ5(~x)µ4(~x) µ2(~x) µ6(~x)µ5(~x) µ6(~x) µ3(~x) , Σ(~x) =σ1(~x) σ4(~x) σ5(~x)σ4(~x) σ2(~x) σ6(~x)σ5(~x) σ6(~x) σ3(~x) , (2.3)where µ l,σ l : Ω→ R; l = 1, ...,6. Under these assumptions, these coefficientsmodel the anisotropic and highly heterogeneous behavior of the medium in thegeophysical problem that is considered. Since in practice these coefficients tendto vary over multiple spatial scales and several orders of magnitude, let us as-sume that the entries of these coefficients are piecewise smooth functions with16jump discontinuities at the material interfaces. These coefficients are referredas the medium parameters throughout this dissertation. The stated PDE systemassumes International System Units and eıωt time dependence.Note that taking the divergence of equation (2.1) leads to∇ · ~B = 0, in Ω, (2.4)which implies that ~B must be divergence free for every point in the domain. Thiscondition is known as the Gauss’s law for magnetic fields [58].As shown in [68, 169], the system (2.1)-(2.2) is typically closed with the naturalboundary conditions given byµ−1(~x)~B(~x)×~n =~0, ∀~x ∈ ∂Ω, (2.5)or with the non-homogeneous Dirichlet boundary conditions given by~E(~x)×~n = ~E0(~x)×~n, ∀~x ∈ ∂Ω, (2.6)were ∂Ω denotes the boundary of Ω,~n denotes the unit outward-pointing normalvector to ∂Ω, and ~E0 specifies the tangential components of ~E at ∂Ω. However,more general boundary conditions can be imposed to the Maxwell’s system (2.1)-(2.2) as discussed in [104, 169]. For the sake of simplicity in the exposition, thediscussion is limited to the boundary conditions (2.5) and (2.6).For realistic geophysical settings with highly heterogeneous medium param-eters, the quasi-static Maxwell’s equations do not have analytical solutions [169].In such cases, a numerical solution is computed. The next section provides anoverview of the mimetic finite volume discretization method, which is the dis-cretization method I use in this dissertation to do so.2.4 The Mimetic Finite Volume DiscretizationIn this section, I provide an overview of the mimetic finite volume discretizationmethod (MFV). Full derivation details can be found in [68, 93–98]. Since theMFV method provides a conservative, consistent and stable discretization of the17Maxwell’s equations introduced in Section 2.3, I use it as a building block to de-velop the multiscale and upscaling approaches proposed in this work. However,edge-based FE or FD discretization methods as the ones proposed in [104, 121,131, 135, 149, 165], which also provide a conservative, consistent and stablediscretization for Maxwell’s equations, can be used as well.The MFV is an extension of Yee’s FD method ([178]) that constructs discretecurl, divergence and gradient operators satisfying discrete analogs of the maintheorems of vector calculus involving such operators. For example,∇× (∇φ)=~0,where φ is a scalar potential, and ∇ · (∇× ~A) = 0, where ~A is a vector potential.Therefore, the discrete differential operators obtained with MFV do not have spu-rious solutions and the divergence-free condition for the magnetic field (2.4) isautomatically satisfied in the discrete setting. In addition, MFV models the correctbehavior of the tangential and normal components of ~E and ~B through material in-terfaces (i.e., the tangential components of ~E and the normal components of ~B arecontinuous through material interfaces, and the normal components of ~E and thetangential components of ~B are discontinuous through material interfaces [169]).Due to these reasons, the discretization is called mimetic as it preserves (mimics)the underlying relevant mathematical and physical properties of the quasi-staticMaxwell’s equations. Furthermore, MFV leads to sparse and symmetric linearsystems of equations and it can be implemented on general polygonal and poly-hedral meshes [98].Following the guidelines provided by Hyman and Shashkov ([93, 94]), the MFVmethod begins by considering the weak form of the Maxwell’s system (2.1)-(2.2):(∇× ~E,~F)+ ıω(~B,~F) = 0, (2.7)(∇× µ−1~B, ~W )− (Σ~E, ~W ) = (~Js, ~W ), (2.8)where ~F ∈H (div;Ω) and ~W ∈H (curl;Ω) are arbitrary test functions;H (div;Ω)and H (curl;Ω) are the Hilbert spaces of square-integrable vector functions onΩ with square-integrable divergence and curl, respectively; and (·, ·) denotesthe inner product given by (~P, ~Q) =∫ΩPxQx + PyQy + PzQz dV . In particular,~E ∈ H (curl;Ω) and ~B ∈ H (div;Ω). For a more detailed description of thearbitrary test functions ~F and ~W used here, the interested reader is referred to18[93, 94, 131].As shown in [68], after integrating by parts the term (∇× µ−1~B, ~W ) in equa-tion (2.8) and applying the natural boundary conditions (2.5), the boundary termvanishes and equation (2.8) becomes(µ−1~B,∇× ~W )− (Σ~E, ~W ) = (~Js, ~W ). (2.9)The resulting weak system formed by equations (2.7) and (2.9) is more convenientto work with. It only requires the tangential components of ~E and ~W be differen-tiable, and eliminates the differentiability condition on the tangential componentsof ~B, which may be discontinuous across material interfaces. Furthermore, thediscretization of the anisotropic discontinuous coefficients, Σ and µ , will be donevia the numerical approximation of their corresponding inner products; as we willsee, this allows an adequate treatment in such cases.Next, the method proceeds by numerically approximating the differential op-erators and the inner products in the weak system formed by equations (2.7) and(2.9). To elaborate further, let us first introduce some mathematical notation. Forsimplicity in the exposition, let us assume Ω to be a cuboid domain that is dis-cretized with a staggered tensor mesh,M. In particular,M=∪nx,ny,nzi, j,k=1Ωhi, j,k, wherenx, ny and nz are the number of cells along the x, y and z axis, respectively; andΩhi, j,k denotes the (i,j,k)th cell. The lengths of Ωhi, j,k along the x, y and z axis arehxi , hyj and hzk, respectively. Let us number the cell center of Ωhi, j,k as (i, j,k), thex, y and z edges of Ωhi, j,k as (i, j± 12 ,k± 12), (i± 12 , j,k± 12) and (i± 12 , j± 12 ,k),respectively; and the x, y and z faces of Ωhi, j,k as (i± 12 , j,k), (i, j± 12 ,k) and(i, j,k± 12), respectively. Let us discretize ~E and ~Js on the edges, ~B on the faces,and the PDE coefficients µ and Σ at the cell centers. Figure 2.1 shows a controlvolume cell with the allocation of these variables. According to [58, 68], edge vari-ables can be physically interpreted as (electric) fields, as they represent a (line)force moving on a path around a surface. Similarly, face variables can be physi-cally interpreted as (magnetic) fluxes, as they represent a (magnetic) flux inwardsor outwards the cell. The physical meaning of such variables explains the choiceof discretization location. The corresponding grid functions onM for the variablesmentioned before are denoted as e, js, b, µ and Σ, respectively.19xzyΩhi, j,kxi− 12 xi+ 12y j+ 12y j− 12zk+ 12zk− 12Bxi+ 12 , j,kBzi, j,k+ 12Byi, j− 12 ,k{Jxs ,Ex}i, j− 12 ,k− 12{Jzs ,Ez}i− 12 , j− 12 ,k{Jyz ,Ey}i+ 12 , j,k− 12Figure 2.1: Control volume cell showing the staggered discretization for ~Eand ~Js on its edges, ~B on its faces, and the medium parameters µ andΣ at the cell centers.To discretize the curl operator on each of the faces of the cell Ωhi, j,k, let us useits integral form as given in [145], which is convenient as it provides a geometrical,coordinate-free expression for the definition of this operator. That is,(∇× ~E) ·~nS = limAS→0(1AS∫S∇× ~E ·dS)= limAS→0(1AS∮∂S~E ·dl), (2.10)where S denotes a face of Ωhi, j,k, ~nS denotes the outward unit normal vector toS, ∂S denotes the boundary of S, and AS denotes the area of S. In the aboveexpression, the Stoke’s theorem was applied ([145]).To compute the integrals in (2.10), let us use the midpoint quadrature rule([6]). For example, the discretization for the (i+ 12 , j,k)th face, which is denotedas Si+ 12 , j,k, can be written as(∇× ~E) ·~nSi+ 12 , j,k≈−hyj(Eyi+ 12 , j,k+12−Eyi+ 12 , j,k− 12)+hzk(Ezi+ 12 , j+12 ,k−Ezi+ 12 , j− 12 ,k)hyjhzk.(2.11)The discretization for the rest of the faces of Ωhi, j,k can be done in a similar way.20The above computation is performed for every cell and its corresponding faces inthe staggered meshM. The result can be expressed in matrix form as follows∇× ~E ≈ CURLe = S−1CLe, (2.12)where S is a diagonal matrix that contains the area of each face in our mesh, Lis a diagonal matrix that contains the length of each edge in our mesh, and C isa matrix that contains the values 0 and ±1 that indicates the mesh connectivity(i.e., it encodes the signs in equation (2.11)).The discretization of the divergence and gradient operators can be done ina similar manner using their corresponding geometrical, coordinate-free definition(as given in [145]). The interested reader is referred to [68, 94, 96] for more detailson how to do so. In such references, the authors demonstrate that the discreteoperators obtained in this way satisfy thatCURLGRAD= 0M, DIVCURL= 0M, (2.13)where GRAD denotes the discrete gradient operator, DIV denotes the discretedivergence operator, and 0M is a matrix with all entries equal to zero. Further-more, they demonstrate that GRAD spans the non-trivial null space of CURL.These discrete relationships are known as the mimetic properties of the discreteoperators.To complete the discretization procedure, we need to approximate the innerproducts (Σ~E, ~W ), (µ−1~B,∇× ~W ), (~Js, ~W ) and (~B,~F). As we will see, the approx-imation of (∇× ~E,~F) will be done in an analogous way as the approximation of(~B,~F). We continue to follow Hyman and Shashkov ([93, 95–97]), and use low-order quadrature formulas to do so. One of the main differences between MFVand an edge-based FE discretization for Maxwell’s equations ([104, 131]) is thatMFV uses low-order quadrature formulas (e.g. midpoint or trapezoidal rules) forthe numerical integration of the corresponding inner products, whereas FE usesGauss numerical integration rules.First, let us consider the approximation of (Σ~E, ~W ). Recall that the x, y and zcomponents of ~E and ~W are discretized on the edges and that Σ is discretized at21the cell center of each cell (Figure 2.1), which implies that Σ is considered to takea constant value within each cell. Since Σ is a SPD matrix of the form (2.3), thecomputation of the integrand Σ~E ·~W involves cross-term products of the tangentialcomponents of ~E and ~W that are not discretized at the same location; thus no low-order quadrature formulas (e.g. midpoint or trapezoidal rules) can be applied inthis case. To overcome this problem, one first projects the tangential componentsof ~E and ~W to the nearest node in the cell; then one computes Σ~E · ~W . Thisprocess is referred to as nearest-neighbor interpolation ([68]) and it is illustratedin Figure 2.2(a). After applying this process, one obtains an approximation to(Σ~E, ~W ) as follows∫ΩΣ~E · ~W dV ≈ w>(188∑n=1Pen>V12M(Σ)V12Pen)e= w>Me(Σ)e, (2.14)where V is a diagonal matrix with the volume of each cell in the mesh, Pen, n =1, . . . ,8, are the matrices that project the three adjacent components of the electricfield to the nth node of each cell, and M(Σ) is given byM(Σ) =diag (σ1) diag (σ 4) diag (σ 5)diag (σ 4) diag (σ 2) diag (σ 6)diag (σ 5) diag (σ 6) diag (σ 3) , (2.15)where each diag (σ i) represents a diagonal matrix containing the grid function ofσ i; i = 1, ...,6, respectively. Let us define Me(Σ) in (2.14) as the edge mass ma-trix. Note that when Σ is a diagonal matrix, the edge mass matrices is identicalto the one obtained by discretizing (Σ~E, ~W ) using low-order quadrature formu-las [70]. The described discretization method naturally extends the use of Yeemeshes (i.e., staggered tensor meshes) to the anisotropic case. Its advantageis that it only requires integration over one cell at a time; the averaging (2.14)is a result of the integration over the domain. This is similar to FE methods,where averaging is performed implicitly through the assembly of the mass matrix[104, 131].Next, let us approximate (µ−1~B,∇× ~W ). Since µ is a SPD matrix of theform (2.3), computing the integrand µ−1~B ·∇× ~W leads to a similar problem as22Figure 2.2: Schematic representation of the nearest-neighbor interpolationprocess for approximating (a) the edge inner product (Σ~E, ~W ) and (b)the face inner product (µ−1~B,~F) in the cuboid cell Ωhi, j,k.the one described before for the numerical integration of (Σ~E, ~W ). In this case,one uses the nearest-neighbor interpolation method on the faces of M and thediscrete expression of the curl operator obtained in (2.12) in order to approximate(µ−1~B,∇× ~W ), which leads to∫Ωµ−1~B · (∇× ~W )dV ≈ w>CURL>(188∑n=1Pfn>V12M(µ−1)V12Pfn)b= w>CURL>Mf(µ−1)b, (2.16)where Pfn; n= 1, ...,8 are the matrices that project the three adjacent componentsof the magnetic flux to the nth node of each cell (see illustration of this process inFigure 2.2(b)), V is defined as before in (2.14), and M(µ−1) is given byM(µ−1) =diag (µ−11) diag (µ−14) diag (µ−15)diag (µ−14) diag (µ−12) diag (µ−16)diag (µ−15) diag (µ−16) diag (µ−13) , (2.17)where each diag (µ−1i) represents a diagonal matrix containing the grid func-tion of the ith entrance of the inverse of the matrix µ given by (2.3), denoted asµ−1i; i = 1, ...,6, respectively. Mf(µ−1) in (2.16) is defined as the face massmatrix.23Next, let us approximate (~Js, ~W ). Since the three tangential components ofthe integrand ~Js ·~W are placed at the same location for every cell of the mesh (i.e.,at the edges), one uses a combination of the trapezoidal and midpoint quadraturerules ([6]) to approximate this inner product, which leads to∫Ω~Js · ~W dV ≈ w>(diag (Acce )>v)js, (2.18)where Acce is a sparse matrix that averages edge variables into the cell centersand it contains only 14 where the edge variables are averaged ([68]), and v is thevector of cell volumes arising from applying the quadrature rules on each cell ofthe mesh.Next, let us approximate (~B,~F). Since the three normal components of theintegrand ~B · ~F are placed at the same location for every cell of the mesh (i.e.,at the faces), one uses the midpoint quadrature rule to approximate this innerproduct, which leads to∫Ω~B ·~F dV ≈ f>(diag (Accf )>v)b, (2.19)where Accf is a sparse matrix that averages face variables into the cell centers andit contains only 12 where the face variables are averaged [68], and v is defined asbefore in (2.18). The approximation of the inner product (∇× ~E,~F) can be donein an analogous way as for (~B,~F) and also uses the discrete expression of the curloperator that was obtained in (2.12). Doing so results in the following expression∫Ω∇× ~E ·~F dV ≈ f>(diag (Accf )>v)CURLe. (2.20)Combining equations (2.14), (2.16), (2.18), (2.19), (2.20) and applying somestandard mathematical manipulations, one obtains the following discrete analogto the quasi-static Maxwell’s equationsCURLe+ ıωb = 0, (2.21)CURL>Mf(µ−1)b−Me(Σ)e = diag((Acce )>v)js =: q, (2.22)where 0 is a column vector with all its entries equal to zero. Now, from equation24(2.21), one can obtain a discrete analog to the magnetic flux as followsb=− 1ıωCURLe. (2.23)Note that multiplying (2.23) times DIV and using (2.13) implies thatDIVb= 0; (2.24)that is, one obtains the discrete analog to Gauss’s law for the magnetic field (2.4).Now, substituting (2.23) into (2.22) one obtains the following system of linearequations in terms of eA(Σ)e=−ıωq, (2.25)whereA(Σ) = CURL>Mf(µ−1)CURL+ ıωMe(Σ). (2.26)The matrix A(Σ) is complex, sparse, symmetric, and, in practice, it tends tobe severely ill-conditioned. Section 2.5 provides an overview of the direct anditerative solvers that can be used to solve this system.As shown in [68], to impose the non-homogeneous Dirichlet boundary con-ditions (2.6), which imply the values of the tangential components of the electricfield at the boundary are known, the matrix A(Σ) and the vectors e and q fromequation (2.25) are reordered into interior edges (ie) and boundary edges (be).Thus, the system to be solved in terms of the unknown eie isAie,ieeie =−(ıωqie+Aie,beebe), (2.27)where Aie,ie, Aie,be, and qie represent the corresponding partitions of the matrixA(Σ) and the vector q of the system (2.25), eie is the discretized electric field atthe interior edges, and ebe is the discretized electric field at the boundary.As shown in [68, 96], the MFV method is second order accurate assuming anorthogonal mesh and that the PDE coefficients are smooth or piecewise constant,which is our case. The works in [124, 131] propose discretization schemes forMaxwell’s equations with higher orders of accuracy. However, these schemes aremore complicated and more computationally expensive to solve. The order of25accuracy provided by MFV is sufficient for the type of large-scale geophysical EMsimulations we are interested in this study.The ideas presented in this section have been extended to OcTree meshes(for details see [72, 86]) and to logically orthogonal meshes (for details see [68]).Let us finish this section by providing a mimetic discrete approximation for theelectric current density ~J, as it will be used later in Chapter 3. The continuousform of Ohm’s law ([169]) states that~J = Σ~E. (2.28)In order to obtain a discrete analog of the electric current density, j, one firstinterpolates the tangential components of ~E to the cell centers of our mesh. Next,one performs a point-wise multiplication of such components with the conductivityto obtain ~J at the cell centers. Finally, one interpolates the three componentsof ~J to its normal directions (i.e., to the cell faces). This procedure results ina conservative interpolation for ~J that preserves its continuous properties in thediscrete setting. The grid function of ~J is denoted as j, which is given byj= AfccMe(Σ)Acce e, (2.29)where Acce represents the averaging matrix taking information from edges to cellcenters as defined in (2.18), Me(Σ) is given by (2.15) and Afcc is a sparse matrixthat averages cell-center variables into the faces and it contains only 12 where thecell-center variables are averaged ([68]).2.5 SolversThis section provides an overview of the direct and iterative solvers used to solvethe large, complex, symmetric, and sparse linear system of algebraic equations(2.25) that result from the MFV discretization of the quasi-static Maxwell’s equa-tions.Before discussing the solver alternatives for the system of equations (2.25),a brief discussion about the conditioning of such system of equations is given.A severe ill-conditioning problem can arise in practice when solving the system26of equations (2.25) for the cases where the EM survey configuration considersvery low frequencies, i.e., ω ≈ 0; and/or when the geophysical setting considersvery low conductivity values (e.g. when the setting considers air, whose electricalconductivity is 0 S/m, but for simulation purposes this conductivity value is typicallyassumed to be 10−8 S/m). For such cases, the matrix system (2.26) is close tobe singular and its condition number is expected to be quite large.The ill-conditioning of the system matrix (2.26) for the cases described abovecan be primarily explained by the fact that the discrete curl operator obtainedusing the MFV discretization has a non-trivial null space, namely the discretegradient operator (see equation (2.13)). This feature is expected from using amimetic method to discretize the quasi-static Maxwell’s system. That is, the dis-crete system is expected to inherit the near-singularity of the underlying continu-ous Maxwell’s system [5, 68].More generally, the conditioning of the system matrix (2.26) can be explainedby analyzing separately its two parts. The work in [40, 68] shows that the eigen-values of the term Me(Σ) in equation (2.26) are bounded by the minimum andmaximum values (assuming volume scaling) of the electrical conductivity consid-ered in the geophysical setting. On the other hand, the largest eigenvalue of theterm CURL>Mf(µ−1)CURL is not bounded. Such term comes from the MFVdiscretization of the differential operator ∇× (µ−1∇×). As one refines the mesh,its largest eigenvalue goes to infinity, while its smallest eigenvalue is zero (due tothe non-trivial null space of the curl operator) (cf. [68]).Once the conditioning of the system (2.25) has been discussed, lets discussthe solver alternatives for solving such system.Iterative solvers are the natural choice to solve the large, sparse and sym-metric system of equations (2.25). They lead to substantial savings on timeand memory usage, as they use optimized (sparse) matrix-vector product mul-tiplications at each iteration and do not require to store the matrix of the system[146]. The literature suggests using BiCGSTAB, MINRES or GMRES solvers forthe system (2.25) [68, 134, 146]. Such solvers can experience extremely slowconvergence rates or non-convergence at all when the system (2.25) is extremelyill-conditioned (e.g. when simulation setting considers very low frequencies andvery low conductivity values) [75]. For such cases, the literature suggests that we27first reformulate the PDE model to obtain a better conditioned system. And sec-ond, that we precondition the system of equations obtained from the reformulationof the problem. Doing so results in faster convergence rates to the solution (e.g.[5, 12, 66, 68, 70, 100, 115, 128, 135, 149, 157, 168, 172]). Below, we brieflydiscuss these two complementary steps:1. The first step is to reformulate the Maxwell’s system (2.1)-(2.2) using theHelmholtz decomposition together with a Coulomb gauge condition suchthat, after applying an adequate discretization method, the resulting ma-trix is better conditioned. The Helmholtz decomposition is a fundamentaltheorem of vector calculus that expresses the electric field as the sum ofa curl-free potential field and a divergence-free potential field [145]. Thisapproach can lead to a saddle point system of equations, which can besolved leveraging its special structure (see [63, 64] and references within).To avoid solving the saddle point system, further algebraic manipulationscan be done to the reformulated system in order to obtain a much betterconditioned system to be solved (see [68, 152] and references within).2. After the issues associated with the null space of the curl have been takencare of by reformulating the PDE system as discussed in step 1, the sec-ond second step is to construct a preconditioner matrix. Such precondi-tioner should transform the reformulated system to be solved into anothersystem with more favorable properties, such as a clustered spectrum, foriterative solution [15, 146]. The literature suggests the use of the follow-ing preconditioners: Jacobi, Symmetric Gauss-Seidel, Incomplete LU andSSOR [68, 100, 134]. As discussed in [68], these preconditioners are easyto apply and work well for a moderate-sized system, if the system is notpoorly conditioned. For much larger ill-conditioned systems, more sophisti-cated preconditioners are needed, such as multigrid [27, 164]. Research onmultigrid solvers and preconditioners for Maxwell’s equations can be foundin [4, 83, 84, 105, 133].Although using iterative methods (complemented with a reformulation of thePDE model as described above) is one of the most efficient ways currently used totackle EM forward simulations, using iterative methods as part of an EM inversion28routine can suffer a major drawback: the system needs to be solved multipletimes for the number of sources and the different frequencies considered in thegeophysical EM survey configuration [68], which can be in the order of thousandsor even millions (e.g. an airborne survey) [144]. In such cases, being able to usedirect solvers to decompose the matrix of the system of equations is the mostcost-effective way to overcome this issue [43, 68]. However, a matrix factorizationrequires us to be able to store the matrix. This requirement is what motivates myinterest on being able to reduce the size of the system of equations, so that directmethods can be used to solve the system. Let us discuss next the options fordirect solvers for the linear system (2.25).Direct solvers are based on a version of Gaussian elimination to compute aLU factorization [43]. Since the factors are stored, these methods require to haveenough CPU memory available. These type of solvers are known to produce morereliable and robust solutions than iterative solvers. The research in [149] pointedout that for settings that use frequencies and conductivities very close to zero, thesystem is severely ill-conditioned and the PDE problem may benefit from beingfirst reformulated as discussed previously.For 1D problems and most 2D problems direct factorization methods are avail-able on most packages and can be used without specialized computing resources.However, for 3D problems it is not always possible to factorize the system withoutspecialized software and hardware. There are some parallel direct solvers pack-ages available, such as MUMPS [3], SuperLU [116], STRUMPACK [62], PARDISO[148], WSMP [65], that can solve problems with millions of unknowns. Thesepackages use optimized memory saving strategies and its matrix factorization isdone in parallel architectures on powerful hardware, large-scale problems can behandled with direct solvers. My research group has used the solver MUMPS forseveral years with satisfactory results, thus the simulation results presented in thisdissertation use such direct solver.2.6 Literature ReviewThis section reviews some of the most successful upscaling and multiscale meth-ods, developed for single-phase flow in porous media problems, that are used29as a building block to develop their counterparts for application to geophysicalEM modeling. In addition, it discusses some of the existing connections betweenmultiscale and multigrid methods.Upscaling and multiscale methods have been extensively studied also forproblems in the fields of materials science, computational mechanics, petroleumengineering and computational water resources [51, 60]. All of these problemsshare the same underlying mathematical model: a Poisson-type equation (i.e., ascalar, linear and elliptic PDE). Since the amount of research for both upscalingand multiscale methods for (linear) elliptic problems is quite extensive, the follow-ing sections only discuss the closest investigations to the upscaling and multiscaleapproaches used in this dissertation. Fewer research along these lines has beendone for the quasi-static Maxwell’s equations discussed in Section 2.3.2.6.1 Upscaling MethodsUpscaling or homogenization methods seek to derive a coarse-scale PDE and topre-compute its coefficients from a given fine-scale PDE model with highly dis-continuous coefficients (see [49, 56, 60, 111, 143, 147, 175] for reviews). Thecoarse-scale PDE coefficients are referred in the literature as upscaled, equiva-lent, averaged, homogenized or effective coefficients [122]. In practice, the setupfor an upscaling method assumes a fine mesh, which accurately discretizes thefine-scale PDE, and a much coarser mesh, which is overlaid on top of the finemesh, where the coarse-scale PDE is ultimately solved. The computation of theupscaled coefficients is done for every coarse-mesh cell. Once the upscaled co-efficients are computed, one can use them to solve the coarse-scale PDE on amuch coarser mesh using fewer computational resources.Using methods from asymptotic homogenization theory ([14]), it was provedthat for a Poisson’s equation that satisfies certain conditions (e.g. scale separationand periodic boundary conditions), the fine- and coarse-scale equations are of thesame form, except that the fine-scale coefficient is replaced by the upscaled one(see [51, 122] for details). In general, by using similar analysis techniques to studyfine-scale heterogeneous structures can result in full tensor upscaled coefficients(see [119, 163] and references within). For example, the investigation in [103]30demonstrates that upscaling two-scale periodic media can result in a full SPDtensor upscaled coefficient, even though the fine-scale coefficient is a scalar.Although technically replacing the fine-scale Poisson’s equation with its anal-ogous coarse-scale equation is only valid under very particular conditions, theliterature provides numerical evidence that justifies such replacement in practice(cf. [49, 122]). In light of this knowledge, when using upscaling methods to solveflow in porous media problems in practice, the focus goes in computing the up-scaled coefficient.The research works in [102, 110, 173, 174] propose similar proofs (under sim-ilar restrictive assumptions) for the case of the quasi-static Maxwell’s equations inthe frequency domain. Based on these studies, in this work, I also assume thatthe fine- and coarse-scale Maxwell’s equations are of the same form and I willfocus on computing the upscaled coefficients of this mathematical model.There is extensive research on how to compute accurate upscaled coefficientsfor flow in porous media problems. Some of the most popular procedures havethe potential to be extended for application to geophysical EM problems, such asanalytical upscaling, flow-based upscaling, and output least-squares upscaling.Analytical UpscalingOne of the simplest ways to compute upscaled coefficients is using an analyti-cal procedure. Analytical procedures seek to derive closed-form expressions forthe upscaled coefficients using averaging principles. It is often the case that anupscaled coefficient is computed by using simple averages of the fine-scale coef-ficient information inside the target coarse-mesh cell. See [49, 56, 60, 141, 175]for reviews of analytical upscaling procedures.Analytical procedures work well when the fine-scale coefficient varies over afixed number of distinctive length scales and it has a particular structure (e.g. alayered medium or a medium with a small correlation length). Similar closed-formexpressions for the upscaled coefficients can be obtained using asymptotic ho-mogenization theory or effective medium theory. Along these two lines, the worksin [28, 102, 103, 110] and in [18, 127, 151, 163] propose analytical proceduresto compute the upscaled electrical conductivity in the context of EM problems,31respectively.Analytical upscaling procedures have the advantage of being simple, verycheap to compute, and accurate in the (limited) cases where the assumptionsare satisfied; however, their underlying fundamental limitation is that these pro-cedures are quite inaccurate for arbitrary fine-scale coefficient variations, whichappear most of the time in practice [49, 56, 122, 132, 175].Flow-based UpscalingA more general and accurate, but more expensive, way to compute upscaledcoefficients for flow in porous media problems is to use a so-called flow-basedprocedure (cf. [45, 49, 170]). This method has its origins in early 1960s. Currently,this is one of the most popular procedures to perform upscaling in the field ofpetroleum engineering.The goal of a flow-based upscaling procedure is to construct an upscaled co-efficient (e.g. permeability) for every coarse-mesh cell by averaging fine-scale flowinformation within the target coarse-mesh cell. Such fine-scale flow is computedby numerically solving a set of local diffusion (steady-state elliptic) boundary-valueproblems (without source term) for a given set of boundary conditions on the tar-get cell. The upscaled coefficient is determined such that the total flow acrossthe coarse cell is preserved as much as possible. These types of procedures canproduce full-tensor upscaled coefficients.Flow-based procedures can be classified as local, extended, or global meth-ods, depending on the size of the computational region used to determine the up-scaled coefficient (see [49, 56, 60, 175] for reviews and discussion). With a localupscaling method, the upscaled permeability is computed using only fine-meshinformation within the target coarse cell. For an extended method, the upscaledpermeability is computed using fine-mesh information within an extended coarse-mesh cell, which includes the target coarse cell and a neighborhood of fine-meshcells around the target coarse cell. Extended methods are more accurate than lo-cal methods mainly because they allow larger heterogeneous connectivities in thetarget cell to be represented more accurately as well as they mitigate the effectof the chosen local boundary conditions, but they are also computationally more32expensive. In a global method, the upscaled permeability for a target coarse cellis computed using the entire fine-mesh model. Global procedures are significantlymore accurate as they can better capture the connectivity effects of the perme-ability throughout the domain, and they do not require local boundary conditions.Global upscaling procedures can be quite computationally expensive for large-scale problems as they require to solve the fine-mesh PDE problem more thanonce. For very large-scale 3D problems, they are still not practical.With a proper setup, flow-based procedures can drastically reduce the costof a fluid flow simulation, while providing an accuracy comparable to the one ob-tained with traditional fine-mesh discretization techniques, such as FE or FV (see[35, 36, 48, 49, 60, 61, 175] for details). However, local and extended flow-basedprocedures can lack accuracy when the local boundary conditions chosen are notproperly chosen, and when heterogeneous connectivities have a spatial extent/s-cale that is larger than the target coarse cell.To overcome the lack of accuracy in local and extended flow-based proce-dures, global procedures can be combined with local or extended flow-basedmethods. Such combined procedures are referred in the literature as local-globalupscaling methods. The investigations in [35–37, 49, 61] are examples of somelocal-global methods that show how global upscaling procedures can be per-formed using a simplified version of the fine-scale PDE in order to obtain a betterset of boundary conditions to setup a local or an extended upscaling procedure.A significant amount of research has been developed in order to improve theflow-based upscaling procedures. Some of the most important research directionsto do so are: a) testing different boundary conditions for the set of local diffusionsproblems that need to be solved for 2D and 3D problems, b) understanding how tobetter capture full-tensor effects derived from the upscaling procedure c) investi-gating how the combination of different upscaling procedures perform for differentproblems, d) proposing flow-based coarse-mesh generation methods and investi-gating the effect of this type of mesh setup in the resulting upscaled coefficients,and e) assessing the quality of the resulting upscaled permeability model. See[35, 36, 49, 50, 60, 61, 112] and references within for details and discussion re-garding these topics.33The Output Least-Squares UpscalingAmong the various existing global upscaling formulations for flow in porous mediaproblems, there is one that can be generalized to solve upscaling problems indifferent contexts: the output least-squares global method (OLS) (for reviews see[35, 49, 56, 175]). The OLS method can be considered as a global procedureas it computes upscaled coefficients by performing computations on the entirefine-mesh model.The OLS method was introduced in early 2000s in the studies by [85, 136],where the authors propose a least-squares formulation for the upscaling problem.All the upscaled permeabilities are computed at once by minimizing the regular-ized least-squares difference between the pressure and the velocity fields gen-erated using the fine and coarse-scale pressure equations. These studies onlyprovided examples for constructing scalar upscaled coefficients for 2D problemswith promising results. However, the OLS method is quite expensive to performas the fine-mesh problem needs to resolved several times.A least-squares formulation can be customized to compute either isotropic orfully anisotropic upscaled coefficients, depending on the definition of the objectivefunction and the purpose of the simulation. Due to the flexibility and generality thatthis formulation offers, I extend and adapt it to solve geophysical EM problems inChapter 3.2.6.2 Multiscale Finite Element and Finite Volume MethodsThere is a considerable body of research proposing multiscale approaches forproblems with features at multiple length scales. Most of the work has been de-veloped for flow in porous media problems. Among the most popular methodsare the FE heterogeneous multiscale method [1, 171], the variational multiscalemethod [92], the generalized FE multiscale method [9–11], and the multiscaleFE/FV methods [51, 89, 101]. The book in [51] provides a comprehensive reviewon this topic.In this dissertation, I focus on the multiscale FE/FV methods develop for linearelliptic problems originally proposed by [74, 89, 101]. This type of multiscale meth-ods have been successful in reducing the size and cost of the computation while34producing accurate solutions similar to that obtained with FE or FV discretizationschemes on a fine mesh.Multiscale FE/FV methods belong to the family of projection-based model re-duction methods, where the fine-mesh system resulting from the discretizationof the PDE model is projected onto a reduced subspace [51, 89]. Projected-based methods approximate the unknown physical responses using a basis ofreduced dimension and project the governing PDE onto a suitably defined low-dimensional subspace [13]. The setup for these method assumes a fine mesh,which accurately discretizes the fine-scale PDE, and a much coarser mesh, whichis overlaid on top of the fine mesh, where the problem is ultimately solved. Multi-scale FE/FV methods construct multiscale basis functions for every coarse-meshcell. The multiscale basis functions are computed by solving sets of local (elliptic)boundary-value problems. This process ensures that the basis captures the fine-scale variations of the PDE coefficients. A global variational formulation couplesthese basis functions to provide an accurate coarse-mesh solution to the problem.According to [51], the multiscale FE method was born in early 1980s fromthe work in [9, 10], where the authors extended the Galerkin FE method for lin-ear steady-state elliptic problems by constructing multiscale basis functions thatdepend on a particular type of multiscale coefficient. By early 1990s, the inves-tigations in [89, 90] generalized the construction of multiscale basis functions forarbitrary coefficients. In particular, these investigations pointed out that the ac-curacy of the method strongly depends on the local boundary conditions chosento construct the multiscale basis functions in the target coarse cell and they pro-posed oversampling techniques to solve the related accuracy issues. By early2000s, multiscale FE methods were extended to non-linear, time-dependent (el-liptic and parabolic) problems and various other global variational formulationsto couple the multiscale basis functions were proposed (see [51] and referencewithin). During this time, the investigations in [79, 81, 101, 160] proposed mul-tiscale FV methods, where a FV global formulation is used and the multiscalebasis functions are computed by discretizing local boundary-value problems ina staggered mesh (as opposed to a nodal mesh as in multiscale FE). By doingso, the multiscale FV method is able to produce mass-conservative, coarse-meshsolutions.35Since early 2010s, some of the most important research directions in the mul-tiscale FE/FV community have been: a) extending the range of applicability ofthese methods to simulate other physical phenomena, b) analyzing error con-vergence and improving the approximation properties of the method for differentsettings, c) generalizing the method by using more basis functions and improvingoversampling techniques, d) developing multilevel multiscale methods as well asdesigning efficient parallel implementations, and e) developing hybrid multiscalemethods using ideas from upscaling techniques (see [39, 50–52, 82, 114] for dis-cussion and details on these topics).Contrary to upscaling methods, which have not been applied to geophysicalEM problems, Haber and Ruthotto ([73]) extended the work by Hou and Wu ([89]),Jenny et al. ([101]), and MacLachlan and Moulton ([123]) that propose multiscaleFE/FV for elliptic equations, to develop a multiscale FV method for the quasi-staticMaxwell’s equations in the frequency domain. However, their work did not includeoversampling and therefore, their method can result in inaccurate solutions due to‘resonance errors’ (i.e., errors that appear when the mesh size and the wavelengthof the small-scale oscillation of the coefficient are similar) [51, 89, 90].In Chapter 4, I extend the oversampling technique proposed by Hou and Wu([89]) developed for elliptic equations, for application to the quasi-static Maxwell’sequations. IIn parallel to the multiscale community, the multigrid community has also beenvery active in proposing solutions to tackle boundary-value problems with multiple-scale features (for reviews see [22, 167] and references within). In the next sec-tion, I discuss some of the connections between the multiscale FE/FV and multi-grid approaches proposed by these two communities.2.6.3 Links between Multiscale FE/FV and Algebraic MultigridMultigrid or multilevel algorithms, introduced in the seminal paper by Brandt in1977 ([21]), provide a methodology to design faster iterative solvers for systemsof algebraic equations. Such a methodology delineates the principles to combinelocal processing on different levels (scales) with inter-level (inter-scale) transferoperators in order to design solvers that can yield linear (optimal) complexity, i.e.,36the computing processing is proportional to the problem size [24]. Since theirinception, this methodology continues to be generalized to tackle a wide rangeof problems. Multigrid can be classified into two main branches: geometric andalgebraic methods (see [24, 164] for details).The design of an efficient algebraic multigrid (AMG) solver requires an effec-tive coarsening process, i.e., the process of generating coarse-level problems (oroperators) and inter-level restriction and interpolation operators. Such operatorsmust satisfy certain weak-approximation properties, i.e., some necessary condi-tions for the iterative process to be convergent [55, 167]. The coarsening processis usually automatic and operator dependent (i.e., it uses information of the matrixof the system to be solved). An effective coarsening process to design robustAMG solvers is the so-called variational Galerkin coarsening [164].Variational Galerkin coarsening automatically generates a hierarchy of coarse-level problems and coarse-to-fine interpolation operators by using minimizationprinciples. For example, for a two-level method, given the finest-level operator(A f , which is the matrix of the system to be solved), a coarse mesh, and an in-terpolation operator (P, which is also computed using variational principles), thecoarse-mesh Galerkin operator (Ac =PT A f P) that minimizes the error in the rangeof the interpolation operator P is readily obtained [109, 164]. In this method, therestriction operator (R) is given by R = PT . The construction of P gives rise to acoarse space defined by the range of P. That is, the construction of interpolationoperators for a hierarchy of levels implies that a sequence of coarse spaces isalso generated [26].On the other hand, (Galerkin) multiscale FE/FV methods, discussed in theprevious section, compute multiscale basis functions by solving local problems oneach coarse-mesh cell. These basis functions are used in the weak formulationof the boundary-value problem of interest in order to create a coarse-mesh dis-cretization operator Ac that has a similar form as the Galerkin operator describedabove (i.e., Ac = PT A f P), where the columns of P contains the assembled mul-tiscale basis functions computed in each coarse-mesh cell. The above outlinesone of the main connections between AMG and multiscale FE/FV methods.37AMG as a Numerical Upscaling ToolThe AMG methodology has also been used as a numerical upscaling tool, i.e., tosolve a coarse-mesh problem (of the form PT A f Pxc = PT b) rather than a largesystem of equations (of the form A f x = b) that results from the FE discretizationof a PDE on a fine mesh (see [167] for a review). Note that in the AMG context,the word ‘upscaling’ has a rather different meaning than the one given in Section2.6.1. Here, it refers to an upscaled model, instead of a model with upscaledcoefficients. Two AMG-based upscaling methods that are related to multiscaleFE/FV methods for elliptic problems are the element-based AMG method and themultilevel multiscale mimetic method.The element-based AMG (AMGe) method, originally proposed to design effi-cient AMG solvers in [26] and later extended for upscaling purposes in [113, 140],uses FE information in order to create a hierarchy of coarse spaces with guar-anteed approximation properties that can be used as FE discretization spaces.The approximation properties of the constructed FE-based coarse spaces refer tothe necessary and sufficient conditions to guarantee proper error estimates whenAMGe is used as a discretization method. For a given level, AMGe computes thecorresponding interpolation operator P by assembling local element-based inter-polation operators that are computed on each algebraically-defined element in thislevel. For each element, the corresponding interpolation operator is computed bysolving a local problem that is formulated using the FE information (including thelocal stiffness matrix) of such element. AMGe focuses on constructing hierarchiesof FE-based coarse spaces for the De Rham complex, i.e., the sequence of H1-conforming, H(curl)-conforming, H(div)-conforming and L2-conforming spaces.The De Rham complex constitutes a convenient sequence of spaces where abroad class of PDEs problems can be handled. The method can work for unstruc-tured (conforming) meshes and high-order FE. AMGe upscaling has been mostlyapplied to elliptic problems (e.g. [38, 114]). To the best of my knowledge, AMGe-upscaling applications for the Maxwell’s equations have not been published yet.The multilevel multiscale mimetic (M3) method proposed in [119, 120] for time-dependent elliptic problems is also inspired in using the AMG methodology as anupscaling tool. M3 uses Galerking coarsening in order to construct a hierarchy of38mimetic FD-based coarse-mesh problems. The accuracy of this coarsening pro-cedure to generate a sequence of coarse-mesh discretizations was first shown in[122, 123], where a variational upscaling method is proposed for the Poisson’sequation. However, this method does not produce locally conservative massfields. In contrast, the M3 is locally mass conservative in all levels and worksfor polyhedral meshes. To the best of my knowlege, extensions of M3 to solveMaxwell’s equations have not been published yet.39Chapter 3An Upscaling Framework for theElectrical ConductivityLook closely at the present you are constructing:it should look like the future you are dreaming. — Alice Walker3.1 OverviewThis chapter1 proposes an upscaling framework for the electrical conductivity ofthe mathematical model introduced in Section 2.3. The goal of this framework isto construct accurate coarse-mesh electrical conductivity models from given fine-mesh ones that can be used for faster simulation of the quasi-static EM responseson a coarse mesh.The main idea behind the proposed framework is to pose upscaling as a pa-rameter estimation problem to be solved on each coarse-mesh cell. As we will seein Section 3.2, this formulation is fundamentally different than others proposed inthe upscaling literature. This chapter starts by introducing the components of theframework. Then, a 1D example that upscales well log conductivity data fromthe Canadian McMurray formation is used to illustrate the general principle be-hind the framework. Afterwards, the upscaling formulation is adapted in order to1This chapter contains extended and revised versions of the material published in [29] and [32].I am the lead author in both of these publications.40obtain a practical method for constructing 3D full-tensor conductivity models. Fi-nally, the performance of the 3D upscaling formulation is demonstrated by usingtwo examples on a single coarse cell, and one synthetic example based on theCanadian Lalor mine. The last example shows the feasibility of combining the pro-posed upscaling method with OcTree meshes in a parallel environment to boostits performance.3.2 A Least-Squares Formulation of the UpscalingProblemThis section introduces the mathematical framework that poses the upscaling pro-cedure as a parameter estimation problem.Remember, the goal is to simulate EM responses for large-scale geophysicalsettings that consider highly heterogeneous media and features varying at mul-tiple spatial scales. Often, the small-scale features have a significant effect onthe measured EM responses (e.g. the thickness of a steel-cased well in Figure1.3). In such cases, the accuracy of the computed EM responses of interest de-pends on our ability to capture the relevant fine-scale features into the simulationmesh. This results on using a large mesh, which leads to solving a huge system ofequations, making the simulation computationally expensive, or even intractable.Here, we propose an upscaling framework for the electrical conductivity thatprovides an alternative to avoid computing the EM responses of interest on avery large and fine mesh. To do so, the framework constructs upscaled electricalconductivities that vary on a coarser spatial scale and that emulate the effect ofthe fine-mesh electrical conductivity in the EM responses of interest. This set ofupscaled electrical conductivities can be used for discretization on a much coarsermesh, thus reducing the size of the system of equations to be solved, and, insome cases, making the simulation doable. The upscaling process is illustratedin Figure 3.1.Note that we cannot simply replace the fine mesh with a coarse mesh in thesimulation procedure. Overlaying a coarse mesh over the fine mesh implies thatseveral fine-mesh cells are captured into each of the coarse-mesh cells. This isillustrated for a single coarse cell in the central part of Figure 3.1. As a result,41Figure 3.1: Upscaling process. Maxwell’s equations are solved over the do-main Ω, where the electrical conductivity can vary over multiple spatialscales and several orders of magnitude. The smallest spatial scaleover which the conductivity varies defines the size of the fine mesh onwhich the model is discretized. The goal of an upscaling procedureis to homogenize the conductivity inside a subregion Ωup of Ω in or-der to construct an upscaled electrical conductivity model suitable forsimulation on a much coarser mesh.a single coarse cell contains a heterogeneous fine-scale conductivity structure.From the continuity conditions of the EM fields across material interfaces [169],this implies that the EM responses within this coarse cell are non-smooth. Inthis case, using a standard discretization technique on the coarse mesh, such asFE or FV, will produce inaccurate approximations to the EM responses becausesuch techniques assume a certain degree of smoothness in the function to bediscretized. To avoid this complication when replacing the fine mesh with a coarsemesh, we first need to homogenize the fine-scale conductivity inside the coarsecells. That is, for each coarse cell, we need to assign a representative quantity forthe heterogeneous, fine-scale conductivity contained in it. Doing so, we ensuresmoothness in the resulting EM responses inside each cell. Afterwards, we can42safely apply traditional discretization techniques on the coarse mesh. We nowproceed to develop the mathematical formulation for the upscaling procedure wepropose.We begin the development of the proposed upscaling framework by consider-ing the quasi-static Maxwell’s equations in the frequency domain subject to nat-ural boundary conditions. These equations were introduced in Section 2.3. Forsimplicity of exposition, we assume that the magnetic permeability (µ) takes itsfree space value, namely µ = µ0 = 4pi×10−7 Vs/Am. We use these quasi-staticMaxwell’s equations to define a fine- and a coarse-scale Maxwell’s problem.The fine-scale problem considers the Maxwell’s equations where the hetero-geneous electrical conductivity varies over small spatial scales. In our case, theterm fine scale refers to the relevant smallest spatial scale over which the con-ductivity varies, such that this scale defines the size of the fine mesh on whichthe model is discretized to obtain an accurate approximation to the EM responses(see left-hand side of Figure 3.1). We use the superscript f to denote the fine-scale Maxwell’s problem given by:∇× ~E f + ıω~B f = ~0, in Ω, (3.1)∇× (µ−10 )~B f −Σ f~E f = ~Js, in Ω, (3.2)(µ−10 )~Bf (~x)×~n =~0, ∀~x ∈ ∂Ω, (3.3)where Ω, ∂Ω, ~E f , ~B f , ~Js, Σ f , ω ,~n and ı are defined as before in Section 2.3.On the other hand, the coarse-scale problem considers the Maxwell’s equa-tions where the electrical conductivity does not vary over small scale hetero-geneities. The research work in [16, 17, 110, 173, 174] also assumes that thecoarse-scale Maxwell’s equation have the same form as the fine-scale Maxwell’sequations except that the coefficients are replaced by upscaled coefficients. Inour case, the term coarse scale refers to the largest spatial scale over which wecan construct upscaled conductivities suitable for accurate simulation on a muchcoarser mesh (see right-hand side of Figure 3.1). We use the superscript c to43denote the coarse-scale Maxwell’s problem given by:∇× ~Ec+ ıω~Bc = ~0, in Ω, (3.4)∇× (µ−10 )~Bc−Σc~Ec = ~Js, in Ω, (3.5)(µ−10 )~Bc(~x)×~n =~0, ∀~x ∈ ∂Ω, (3.6)where Ω, ∂Ω, ~Ec, ~Bc, ~Js, ω ,~n and ı are defined as before in Section 2.3.The fine- and coarse-scale problems use the same source term ~Js as it is inde-pendent of the fine and coarse-scale variation of the conductivity. The upscalingproblem consists on constructing a coarse-scale electrical conductivity Σc (withoutsmall scale heterogeneities) such that the solution of the coarse-scale problem is,in some sense, close to the solution of the fine-scale problem [56, 136].Now, we define the coarse-scale electrical conductivity Σc asΣc(~x;σup) ={σup, if ~x ∈ΩupΣ f (~x), otherwise(3.7)whereΩup⊂Ω is an upscaling region where we aim to homogenize the fine-scaleconductivity (Figure 3.1), and σup is an upscaled electrical conductivity that aimsto capture the effect of the fine-scale heterogeneous conductivity inside Ωup onthe EM responses. In this framework, the definition of σup depends on the contextof a simulation. For example, σup may be given by a positive scalar, a real or evena complex matrix, depending on the complexity and purpose of the simulation.We explore this idea further through examples in Sections 3.3 and 3.5.We continue by rewriting the fine-scale problem (3.1)-(3.2) and the coarse-scale problem (3.4)-(3.5) asL (Σ f (~x))~u f =~q f and L (Σc(~x;σup))~uc =~qc, (3.8)respectively. Here, L represents the Maxwell operator (in matrix form), ~u f =(~B f ,~E f )> and~uc = (~Bc,~Ec)> represent the fine- and coarse-scale EM responses,respectively; and ~q f and ~qc represent the corresponding source and boundaryconditions as given in (3.3) and (3.6), respectively. The fine- and coarse-scale EM44responses are obtained by inverting the Maxwell operator; that is,~u f (Σ f (~x)) =L −1(Σ f (~x))~q f and ~uc(Σc(~x;σup)) =L −1(Σc(~x;σup))~qc. (3.9)In practice, a geophysical EM forward modeling simulation is used to computepredicted data at the receiver locations or over the entire domain. The computa-tion of the predicted EM data can be expressed as the action of a linear functional,P , on the fine- and coarse-scale EM responses as follows~D f (Σ f (~x)) =P~u f (Σ f (~x)) and ~Dc(Σc(~x;σup)) =P~uc(Σc(~x;σup)). (3.10)Throughout this work, we refer to ~D f and ~Dc as the fine- and coarse-scale EMdata, respectively. Note that when P equals the identity operator, (3.10) returnsthe EM responses in the entire simulation domain.To complete the construction of Σc, we need to define a criterion for choosingthe ‘best’ upscaled conductivity σup in the upscaling region Ωup. As mentionedbefore, the coarse-scale conductivity (Σc) should be able to emulate the effect ofthe fine-mesh electrical conductivity in the EM responses of interest. That is, werequire a criterion able to construct a σup such that the coarse- and fine-scaledata are similar. We therefore propose the following definition.Definition 1. Let Σ f be a given fine-scale electrical conductivity model and let~u f = (~B f ,~E f )> be their induced fine-scale EM responses (i.e., fine-scale elec-tric field and magnetic flux) given by (3.9), for a given angular frequency ω andsources including boundary conditions ~q f and ~qc. Let ~D f and ~Dc, be some pre-dicted fine- and coarse-scale EM data as defined in (3.10), then the upscaledelectrical conductivity in the upscaling region Ωup, denoted as σ∗up, is defined asthe solution of the following parameter estimation problem:σ∗up = argminσupc(σup) =12∥∥∥~Dc(Σc(~x;σup))−~D f (Σ f (~x))∥∥∥22. (3.11)We refer to c(σup) as the upscaling criterion.The parameter estimation problem (3.11) satisfies that the number of param-eters is less than (or equal) to the number of data and it is a well-posed least-45squares problem, hence no additional regularization is required. Our formulationis a type of global upscaling procedure within the classification of upscaling tech-niques proposed in [49, 56] for fluid flow problems in porous media. In a globalupscaling procedure, the entire fine-scale model is simulated for the calculation ofthe coarse-scale PDE coefficient(s).The investigations in [85, 136] propose another least-squares upscaling for-mulation for fluid flow problems in porous media. In such investigations, all theupscaled permeabilities are computed at once by minimizing the regularized least-squares difference between the pressure and the velocity fields generated by thefine- and coarse-scale pressure equations.We solve the parameter estimation problem (3.11) numerically. We use thediscretize-then-optimize approach, where in the first stage we discretize the for-ward problem and in the second stage we solve the finite-dimensional discreteoptimization problem. This approach is typically used in the literature to solvePDE-constrained optimization problems in EM Geophysics (see [20, 67, 68], andreferences within). The main advantages of this approach are that it allows foran efficient integration of convex optimization algorithms with advanced PDE dis-cretization techniques and solvers, and that it leads to more consistent solutions.In order to use this approach, we require both a stable discretization and an ap-propriate convex optimization method to solve the discrete version of (3.11). Inthe examples presented in Sections 3.3 and 3.5 we specify the choices we madein each case.Remarks1. The upscaling definition 1 we propose may look rather involved at first; how-ever, an analogy can be drawn from the computation of an apparent con-ductivity (or resistivity) in a Direct Current Resistivity experiment. Although itwould not be used for the purposes of simulation, the apparent conductivitycan be considered as an upscaled quantity. For a given electrode geometry,the apparent conductivity is the homogeneous halfspace conductivity thatproduces an upscaled response to the one observed [169]. Within this up-scaling context, the apparent conductivity corresponds to the quantity σ∗up,46Ωup is the homogeneous earth (in this case, Ωup = Ω), ~q are the sources,~u are the potentials,P projects the fields onto the receiving electrode loca-tions, and the data ~D f and ~Dc are the measured voltages.2. One aspect of the upscaling definition 1 is that the upscaling regionΩup caninclude its surrounding conductivity structure. That is, the upscaled conduc-tivity can correspond to an on-site sampling of the electrical conductivity.We demonstrate that constructing upscaled conductivities in this manneryields to more accurate approximations to their corresponding induced EMresponses in the examples included in Sections 3.3 and 3.5.3. The upscaling formulation 3.11 considers the connection between the par-ticular EM experiment configuration used in the simulation and the type ofEM responses that it produces into the calculation of the upscaled conduc-tivities. This is accomplished both through the source term and boundaryconditions (~q), and the choice of data of interest (~Dc, ~D f ). For instance,the data of interest can be chosen among the electric or magnetic fields orfluxes (i.e., ~E, ~H,~B or ~J), or some combination of them.4. The upscaling formulation 3.11 is very flexible and provides a user-defined,application-specific framework to upscale the electrical conductivity of ourMaxwell’s system. The upscaled conductivity can be a real or a complex-valued scalar, or a real or complex-valued matrix, depending on the com-plexity of the setting, the purpose of the coarse-scale simulation, and the re-quired accuracy of the solution. Furthermore, the upscaling criterion (3.11)need not be based on least-squares. The examples presented in Sections3.3 and 3.5 are designed to demonstrate these features.5. Having an application-specific framework is important because it accountsfor the fact that there is no unique upscaled conductivity suitable for allsimulation purposes. Indeed, constructing a different upscaling criterionby changing the data simulated or sources used to excite the system typi-cally leads to a different upscaled conductivity. We demonstrate this pointin the examples presented in Section 3.5.1. Additionally, the experiments47presented in Sections 3.3 and 3.5.2 demonstrate that by changing the fre-quency of the survey, it is possible to obtain different upscaled conductivityquantities. This can be explained by the fact that different EM survey config-urations have different sensitivity functions and sample the earth differently.Thus, the effect of heterogeneous conductivities in the EM data we measureoften differs from experiment to experiment. For example, this effect can beobserved when considering Direct Current Resistivity surveys, where anapparent conductivity computed from a pole-dipole survey may be differentthan for a dipole-pole survey [144].In the next sections we discuss the upscaling procedure for 1D and 3D EMgeophysical forward modeling and provide examples that demonstrate its perfor-mance.3.3 Numerical Results in 1DWe show the performance of the proposed global upscaling framework on a 1Dexample. The purpose of this example is twofold. First, we show how to upscalea well log electrical conductivity model for both a single frequency and a multi-frequency airborne loop-loop survey using our framework. Second, we demon-strate that when coarse-scale conductivity models consider the survey configura-tion in their construction, they lead to more accurate simulation results.3.3.1 Simulations Using a Single FrequencyWell logging is the process of recording various physical, chemical, electrical, orother properties of the rock/fluid mixtures penetrated by drilling a borehole intothe Earth’s crust [144]. Well logs are recorded in nearly all oil and gas wells andin many mineral and geothermal exploration and development wells. Well logconductivity measurements2 have high resolution, with samples every few cen-timeters. However, if these conductivity measurements are to be used for earth-scale simulations or inversions, coarse-scale conductivity models, defined on theorder of meters, are needed. This context provides the perfect scenario to test2In reality, well log data considers resistivity measurements. Since electrical conductivity is thereciprocal value of resistivity, we can obtain conductivity measurements using this relationship.48our upscaling framework and to offer a potential alternative to upscale well logconductivity data in practice. Most of the times in practice, log measurements aresimply averaged to obtain effective petrophysical properties, but as we see next,this may lead to large errors.We use a fine-scale electrical conductivity model, Σ f , given by an induc-tion resistivity log from the McMurray/Wabiskaw Oil Sands deposit well log publicdatabase [176]. The McMurray formation is located in Northern Alberta, Canada,and the log used is shown in Figure 3.2(a). Observe that the interval over whichconductivity samples were taken is 80 m, thus we take this as the simulationdomain, i.e., Ω = [0,80] m. In addition, observe that the electrical conductivityranges over four orders of magnitude. The log chosen has 320 measurementstotal, with a measurement taken every 25 cm. Hence, we defined a uniform finemesh whose thickness is consistent with this scale.We consider a standard airborne survey configuration ([144]), with a frequencyof 300 Hz and a horizontal coplanar arrangement for the source receiver-pair. Thesource-receiver pair is located at a height of 40 m above the earth’s surface andhas a separation of 8.1 m. The source produces a magnetic field, which inducescurrents in the earth, producing secondary magnetic fields, which we measure atthe receiver.We aim to simulate the magnitude of the magnetic field (H-field) for the givensurvey configuration using a coarse-scale conductivity model that can be dis-cretized using a much coarser mesh. From now on, we refer to the coarse-scaleconductivity model as the coarse-mesh conductivity model.To construct a coarse-mesh conductivity model that varies on the meter scaleusing the proposed upscaling framework, we need to choose: (a) a suitablecoarse mesh, (b) the type of upscaled quantity to be constructed, and (c) anupscaling criterion.For the coarse mesh, we consider a uniform mesh nested in the fine meshwith 10 m thickness for each coarse layer. Hence, inside the simulation domain wehave eight coarse layers, where each of them contains a large range of fine-scaleconductivity variation to be upscaled, see Figure 3.2(a). We base the estimatefor the proper cell sizes of the coarse mesh on the skin depth ([169]). Practi-cal experience on mesh design for EM problems reported in [68] suggests that49the smallest cell size in the mesh should be a quarter the minimum skin depth.We consider the median conductivity value (2.3×10−2 S/m) and calculate a skindepths of 191.5 m for the frequency of 300 Hz. Therefore, using cell sizes of 10m is sufficient to capture the decaying behavior of the EM fields in this setup.We assume the upscaled conductivity inside each coarse layer is given bya real positive scalar. This is a common assumption in practice when the for-ward modeling code only handles isotropic conductivity models. A more gen-eral assumption is to consider the upscaled conductivity to be a matrix; how-ever, this would require the forward modeling code to be capable of incorporatinganisotropy, which is not always the case. We will elaborate further on 3D caseswhich incorporate anisotropy in Sections 3.4 and 3.5.Since we aim to simulate the magnitude of the magnetic field at the receiver lo-cation, we take this as the EM data to be matched in the upscaling criterion (3.11).By doing so, we have specified the source term (including boundary conditions)and the data of interest, which define the necessary elements of the upscalingcriterion for this example.Once the necessary elements to setup the framework are selected, the con-struction of the coarse-mesh conductivity model is completed by solving a pa-rameter estimation problem for every coarse layer separately. The set of eightupscaled conductivities will form the coarse-mesh conductivity model.As mentioned in Section 3.2, we use the discretize-then-optimize approachto solve each of the eight parameter estimation problem (3.11). To compute theupscaled conductivity on a given coarse layer, we use the EM1DFMfwd code de-veloped by members of the Geophysical Inversion Facility to forward model thenecessary fine and coarse H-field data. The EM1DFMfwd code implements thematrix propagation approach proposed in [57]. Since we have a single H-fielddatum and we are inverting for a single scalar, our parameter estimation problemis well-defined for this example. To solve the discrete version of the optimizationproblem (3.11), we use the MATLAB function fminbnd [162]. Such function imple-ments a standard derivative-free minimization method for single-variable functionson a fixed interval (for more details see [25]). We use a derivative-free optimiza-tion method as an explicit calculation of the derivatives of the upscaling criterionare not available in this case. The resulting coarse-mesh conductivity model is50(a) (b)Figure 3.2: Induction resistivity log: AA-05-01-096-11w4-0 from the McMur-ray/Wabiskaw Oil Sands deposit well log database. (a) Discrete fine-scale conductivity model. Each diamond represents a conductivityvalue on a uniform fine mesh of thickness 25 cm. Each straight line rep-resents a coarse layer of a uniform mesh of thickness 10 m. The setupconsiders 320 fine layers and 8 coarse layers. (b) Resulting coarse-mesh conductivity models after applying four upscaling procedures: 1Dnumerical upscaling procedure (red solid line), and arithmetic (greendot line), geometric (magenta dot dash line), and harmonic (blue dashline) averages.plotted with a red solid line in Figure 3.2(b).Note that calculating an upscaled conductivity following the procedure outlinedabove requires using the fine-mesh conductivity model, which corresponds in theupscaling literature for flow applications to perform a global upscaling procedure[49, 56]. This is a serious drawback to perform 3D simulations in practice. In51Section 3.4 we show how to modify our upscaling framework for 3D settings toavoid this complication.To compare our method to other 1D average-based upscaling methods, weconstruct coarse-mesh conductivity models using volume-arithmetic, -geometricand -harmonic averaging of the fine-mesh conductivity inside each coarse layer.The resulting coarse-mesh conductivity models are shown in Figure 3.2(b). Fromthis figure, we observe that the coarse-mesh conductivity model produced by ourproposed upscaling procedure resembles the coarse-mesh model produced byvolume-arithmetic averaging for most of the coarse cells.To judge the quality of the various coarse-mesh conductivity models as com-pared to the fine-mesh conductivity model, we use each of the models to forwardmodel a H-field datum at the receiver location using the EM1DFMfwd code onthe coarse and fine mesh, respectively. Table 3.1 shows the magnitude of the re-sulting H-field datum for each of the conductivity models and their correspondingrelative errors. The relative error is computed as the ratio of the absolute value ofthe difference in magnitude of the fine- and coarse-mesh datum to the absolutevalue of the fine-mesh datum in magnitude.The results in Table 3.1 demonstrate that the proposed upscaling formulationconstructed an optimal coarse-mesh conductivity model, in the sense of equation(3.11), for the airborne survey configuration given and the surrounding conduc-tivity structures, as it yields the smallest relative error in the approximation of theH-field datum of interest. That is, our method gives and optimal on-site predictionof the upscaled conductivity, as defined in equation (3.11), that is fundamentallydifferent from the one given by the other average-based homogenization methodspresented. Note that by upscaling the conductivity model we can reduce the meshsize from using 320 cells to use only 8 cells without sacrificing much accuracy.In the next section, we demonstrate the effect of considering multiple frequen-cies in the construction of the coarse-mesh conductivity model.3.3.2 Simulations Using Multiple FrequenciesIn this case, we want to construct coarse-mesh conductivity models to simulateH-field measurements at five frequencies logarithmically equispaced in the range52Table 3.1: Magnitude of the magnetic field (H-field) datum and relative errorsresulting from forward modeling using the fine-mesh and four coarse-mesh conductivity models. Note: %∗ denotes percentage of primaryfield.ConductivitymodelMagnitude ofH-field (%∗)Relative error(percent)No. of layersin meshFine 0.0549 —– 320Upscaled 0.0536 2.37 8Arithmetic 0.0468 14.74 8Geometric 0.0102 81.42 8Harmonic 0.0043 92.22 8from 10 to 30,000 Hz (i.e., we use the frequencies of 10, 74, 547, 4,053 and30,000 Hz). To do so, we use the same setup for the fine and coarse meshes,and the same airborne-style survey configuration as described in the previoussection.Considering the median conductivity value (2.3×10−2 S/m), the skin depthsfor the frequencies of 10, 74, 547, 4,053 and 30,000 Hz are roughly 1,049, 386,142, 52 and 19 m, respectively. Thus, using cells sizes of 10 m continues to besufficient to capture the behavior of the magnetic fields in this setup for the firstfour frequencies. We challenge the upscaling method with a larger cell size forthe last frequency.We applied the proposed upscaling procedure by optimizing a discrete ver-sion of (3.11) for each of the eight coarse layers separately, and each individ-ual frequency, as described in the previous section. The resulting coarse-meshconductivity models are shown in Figure 3.3(a). Note that for each frequency,we obtained a different coarse conductivity model. This follows from the factthat the upscaled conductivity model is tailored to match the magnetic field de-termined by the survey parameters. These parameters influence the sensitivityof the magnetic field to the conductivity structure. Hence, varying any of theseparameters alters how the conductivity structure is sampled. As a result, theupscaled conductivity model may take on different values depending on the ex-perimental setting, demonstrating that the proposed upscaling approach provides53a user-defined, application-specific framework. These results also imply that theupscaled conductivity changes as a function of frequency, demonstrating that fre-quency dependence on the coarse-scale considered may arise as a result of local,fine-scale heterogeneity.To conclude this example, we construct coarse-mesh conductivity models us-ing volume-arithmetic, -geometric and -harmonic averaging of the fine-mesh con-ductivity inside each coarse layer. For each of these coarse conductivity models,H-field data were then simulated on the coarse mesh using the given airbornesurvey configuration using the EM1DFMfwd code. We compare the resulting H-field data with those computed using the fine-mesh conductivity model for eachfrequency in Figure 3.3(b). The H-field data shown are given in percentage of themagnitude of the primary field. Table 3.2 shows the relative errors obtained foreach case. The relative error is computed as the ratio of the Euclidean norm ofthe difference in magnitude of the fine- and coarse-mesh data to the Euclideannorm of the fine-mesh data in magnitude. From this table, we observe that, onceagain, the proposed upscaling approach leads to better approximations to the de-sired H-field data than using the coarse models constructed by average-basedupscaling procedures for all the frequencies used.Table 3.2: Relative errors for four coarse-mesh conductivity models and forfive frequencies.Relative errors for coarse-mesh conductivity modelsFrequency(Hz)Upscaled(percent)Arithmetic(percent)Geometric(percent)Harmonic(percent)10 0.64 9.93 83.46 93.3474 2.78 12.76 83.36 93.23547 8.21 15.55 80.85 92.054,053 11.76 15.80 66.83 85.0930,000 0.84 9.59 32.55 55.30Performing 1D upscaling of well log conductivity measurements, as shown inthis example, provides a practical technique to construct accurate coarse-scaleconductivity models to be used on large-scale 3D simulations or inversions. Ourupscaling procedure resulted in more accurate approximations to the EM responses54(a) (b)Figure 3.3: (a) Coarse-mesh electrical conductivity models obtained by us-ing the proposed upscaling method for different frequencies. The setupconsiders 320 fine layers and 8 coarse layers. (b) Magnitude of mag-netic field for each frequency, in % of primary field (%∗), resulting fromforward modeling using the fine-mesh electrical conductivity model(black solid line), the different coarse-mesh conductivity models dis-played in (a) (blue dash line), and the coarse-mesh models producedby using arithmetic (red plus dot line), geometric (gray circle dot line),and harmonic averages (green square dot line).of interest at the cost of more computational performance, as compared to tradi-tional average-based upscaling procedures. However, for 1D problems, using our(global) upscaling framework may be considered to be more appealing as 1D sim-ulations present minimal computational bottlenecks.The example presented illustrates the general principle behind our upscalingframework by showing its performance in 1D well log data. However, applying55the current upscaling formulation (3.11) to a general 3D setting is not practical. Itrequires simulating the data to be matched (3.10) on the entire fine mesh, whichcan be computationally demanding. In the next section, we address the challengeof creating a practical upscaling approach for a 3D setting.3.4 A Local Upscaling Framework for 3D SimulationsApplying the global upscaling framework introduced in Section 3.2 as shown inSection 3.3 is impractical for 3D simulations. It requires simulations using theentire fine-mesh conductivity model. In this section, we adapt the proposed up-scaling framework for practical application to 3D settings.To create a pragmatic upscaling method, we combine our least-squares for-mulation with the methodology proposed by Durlofsky in [48, 49] for the field ofsimulating flow in heterogeneous porous media to EM modeling. That is, usingsome specialized boundary conditions, we apply the upscaling procedure locallyinstead of globally. In our case, this means that for each coarse-mesh cell, welocally solve a parameter estimation problem to construct an upscaled conductiv-ity. Doing so cell by cell, potentially in parallel, yields the desired coarse-meshconductivity model. This approach enables us to solve several smaller problemsrather than a single large one, making this procedure suitable for tackling large-scale EM problems.We now discuss, in detail, how to locally apply the upscaling framework. Weassume that a given fine-scale conductivity model is discretized at the cell-centersof a 3D fine staggered mesh,Mh. The fine mesh sufficiently captures the relevantconductivity variations in the model. We denote the discrete fine-mesh conduc-tivity model as Σh. We aim to construct a coarse-mesh conductivity model, ΣH ,that is also discretized at the cell-centers of a user-chosen 3D coarse staggeredmesh, MH . Typically, MH is much coarser than Mh. Throughout this section,the superscripts h and H denote dependency on the fine and coarse meshes, re-spectively. The fine and coarse meshes are a union of n fine cells and N coarsecells, respectively. That is,Mh = ∪ni=1Ωhi andMH = ∪Nk=1ΩHk , where N n. Forsimplicity, we also assume that the meshes are nested, that is MH ⊂Mh; how-ever, the argument presented here can be extended to include more general mesh56Figure 3.4: Local upscaling procedure for 3D settings. Left: fine-mesh elec-trical conductivity model and example of nested meshes setup. Right:extended domain (ΩH,extk ) for a given coarse-mesh cell (ΩHk ) and re-sulting anisotropic upscaled electrical conductivity (ΣHk ).setups. A sketch of the mesh setup described is shown in Figure 3.4.Since the goal is to apply the upscaling procedure on each coarse-mesh cellindependently and locally, we need to identify: (a) the upscaling region, (b) thetype of upscaling quantity to be constructed, and (c) the data to be matched ina local version of the parameter estimation problem (3.11). We discuss each ofthese choices for a single coarse cell below.Let us consider a single coarse-mesh cell, ΩHk . The upscaling region cor-responds to ΩHk , which is composed of the fine cells and the fine conductivitystructure it encloses. To construct an upscaled conductivity in ΩHk that takes intoaccount the surrounding conductivity structure (i.e., to preserve the on-site sam-pling feature of the fine-scale conductivity), we embedded ΩHk in an extendeddomain, ΩH,extk . This extended domain ΩH,extk includes ΩHk and a neighborhoodof fine cells (and their corresponding conductivity values) around ΩHk . This is illus-trated in Figure 3.4. The size of the extended domain is user-chosen. We explorethe effect of different extensions in the examples presented in Section 3.5.To better represent most of the existing heterogeneity inside ΩHk , we assumethe upscaled conductivity to be constructed for this cell to be a full SPD matrix(i.e., fully anisotropic). Full-tensor effects generally arise at the coarse scale, even57though the fine-scale conductivity is isotropic [35]. We denote the upscaled con-ductivity in ΩHk as ΣHk . That is, ΣHk is a real SPD matrix that can be parametrizedusing six scalars, ~σk = (σ k1 ,σk2 ,σk3 ,σk4 ,σk5 ,σk6 ) ∈ R6, as follows:ΣHk (~σk) =σk1 σk4 σk5σ k4 σk2 σk6σ k5 σk6 σk3 . (3.12)According to Definition 1 given in Section 3.2, in order to construct ΣHk bysolving a parameter estimation problem, we require some fine- and coarse-scaledata to be matched in the upscaling criterion (3.11). To generate such data, theMaxwell system should be excited by either a source or some boundary condi-tions, see (3.10).Since we apply the upscaling procedure locally, we assume that sources donot reside inside ΩH,extk . Therefore, rather than choosing some local source(s)to induce EM responses, we assume that the system is excited by some non-homogeneous boundary conditions. Such boundary conditions should reflect thebehavior of the EM responses of interest across the boundary of ΩHk , which isdenoted as ∂ΩHk .In principle, the correct boundary conditions can be obtained numerically bysolving the fine-mesh problem; however, they are impractical to compute. Oneremedy for this problem, suggested in [49] (for fluid flow problems), is to use aset of linearly independent boundary conditions. Note that using linear boundaryconditions in the context of the EM problem can be appropriate to model the actionof distant sources on ∂ΩHk as such action can be perceived as a ‘plane wave’.Now, to generate the data, we first generate fine- and coarse-scale EM re-sponses by locally exciting the Maxwell system using a set of twelve linearly in-dependent, non-homogeneous Dirichlet boundary conditions (one per edge ofΩH,extk ). This yields the following set of twelve local problems:∇× ~Ekl + ıω~Bkl = ~0, in ΩH,extk ; (3.13)∇× (µ−10 ~Bkl )−Σk~Ekl = ~0, in ΩH,extk ; (3.14)~Ekl (~x)×~n = ~Φl(~x)×~n, ∀~x ∈ ∂ΩH,extk , (3.15)58with l = 1, . . . ,12. Here, ~x = (x1,x2,x3) ∈ R3, ∂ΩH,extk denotes the boundary ofΩH,extk , and ~Φl is a vector function that denotes the lth Dirichlet boundary conditionas defined in Table 3.3. From Figure 3.5, we observe that each ~Φl takes the value1 along the tangential direction to the lth edge of ΩHk and decays linearly to 0 inthe normal directions to the same edge.Table 3.3: Analytical expressions for the set of linearly independent bound-ary conditions used to locally generate data in the cuboid cell ΩH,extk(one per edge). Note that~x = (x1,x2,x3) ∈ R3.~Φ1(~x) = [x2x3,0,0]> ~Φ5(~x) = [0,x1x3,0]> ~Φ9(~x) = [0,0,x1x2]>~Φ2(~x) = [x3(1− x2),0,0]> ~Φ6(~x) = [0,x3(1− x1),0]> ~Φ10(~x) = [0,0,x2(1− x1)]>~Φ3(~x) = [x2(1− x3),0,0]> ~Φ7(~x) = [0,x1(1− x3),0]> ~Φ11(~x) = [0,0,x1(1− x2)]>~Φ4(~x) = [(1− x2)(1− x3),0,0]> ~Φ8(~x) = [0,(1− x3)(1− x1),0]> ~Φ12(~x) = [0,0,(1− x2)(1− x1)]>The set of boundary conditions used ({~Φl}12l=1) form the natural basis func-tions for linear edge degrees of freedom of a hexahedral finite element [104, 131];they can be used to model general linearly varying EM responses. That is, thisset of boundary conditions allows us to generate a set of linearly independent (lo-cal) data, which can be used to formulate a full-rank, overdetermined parameterestimation problem. Similar choices were proposed in [49, 51] for the problem ofsimulating fluid flow in porous media, where the PDE model is the Poisson equa-tion. Different studies for flow in porous media have shown that different choicesof boundary conditions lead to different upscaled quantities [48, 49, 56, 175].To compute the fine-scale EM responses for the kth local problem, we usethe fine-mesh conductivity contained in ΩH,extk , denoted as Σhk . To compute thecoarse-scale EM responses for the kth local problem, we assign ΣHk as the con-ductivity value of every fine-mesh cell inside ΩH,extk . Using the correspondingconductivity values of every fine cell inside ΩH,extk , we can now forward modeleach of the twelve problems (3.13)-(3.15) to obtain a set of fine and coarse EMresponses, respectively. To obtain each set, we discretize each of these problemsusing the MFV method as outlined in Section 2.4. Other traditional edge-baseddiscretization methods, such as an edge-based FE method [104, 131], can beused for this part as well as long as they provide a consistent, conservative andstable discretization for our Maxwell system. After the discretization, we obtain a59set of twelve discrete electric fields, Ek = {ek1 , . . . ,ek12}, and a set of twelve dis-crete magnetic fluxes, Bk = {bk1 , . . . ,bk12}. Note that each ekl is a vector whoselength equals the number of fine-mesh edges in ΩH,extk , and that each bkl is avector whose length equals the number of fine-mesh faces in ΩH,extk .Figure 3.5: Non-zero component of each vectorial basis function ~Φl (definedin Table 3.3) plotted on a unitary cube.Next, we need to choose which fine and coarse data are to be matched andhow the local version of the upscaling criterion (3.11) will be formulated. Wenote that the heterogeneous fine-scale conductivity structure insideΩHk generatesnon-smooth responses on ∂ΩHk . Since the role of the upscaled conductivity is toemulate the effect of the fine-scale conductivity inside ΩHk on the EM responsesat ∂ΩHk , for us it makes sense to consider the data to be either the integral of theelectric field over the twelve edges of ∂ΩHk , or the integral of the magnetic flux60over the six faces of ∂ΩHk . We refer to these data as the total electric fields ortotal magnetic fluxes, respectively.To be more specific, we define the total electric field data asdlm =∫edgem~Ekl ·~τedgem d`; l = 1, . . . ,12, m = 1, . . . ,12, (3.16a)and the total magnetic flux data asdl j =∫face j~Bkl ·~nface j dS; l = 1, . . . ,12, j = 1, . . . ,6, (3.16b)where edgem represents the mth edge of ∂ΩHk , ~τedgem denotes the unit tangentvector to edgem, face j represents the jth face of ∂ΩHk , and ~nface j represents theunit outward-pointing normal vector to face j. The data can be calculated by nu-merically integrating (3.16a) and (3.16b) using the set of discrete fields Ek orfluxes Bk, respectively. Remember that Ek and Bk are discretized on the fine-mesh edges and fine-mesh faces inside ΩH,extk , respectively.As before, the choice of the data to be matched during the local upscaling pro-cedure depends on the context of a given simulation and it significantly influencesthe construction of the upscaled conductivity. We will demonstrate both of thisaspects in the examples presented in Section 3.5. These phenomena have beenalso observed on upscaling procedures for fluid flow problems in porous media[47, 49, 56, 60, 85, 136, 175].Finally, we formulate the local version of the discrete parameter estimationproblem (3.11) to be solved in order to construct ΣHk as follows:~σoptk = arg min~σk∈R6c(~σk) =1212∑l=1∥∥∥dl(ΣHk (~σk))−dl(Σhk)∥∥∥22subject to ΣHk (~σk) is SPD.(3.17)Here, dl(Σhk) denotes the 12×1 vector whose entries are the lth total electric fielddata (3.16a) or the 6×1 vector whose entries are the lth total magnetic flux data(3.16b) computed for the fine-mesh conductivity in ΩH,extk (that is, Σhk), respec-tively. Analogously, dl(ΣHk (~σk)) denotes the vector whose entries are the lth total61electric field data (3.16a) or total magnetic flux data (3.16b) computed for the up-scaled conductivityΣHk (~σk), given by (3.12), extended toΩH,extk , respectively. Thatis, the anisotropic upscaled conductivity for the cellΩHk is given byΣHk =ΣHk (~σoptk ),as parametrized in (3.12). We refer to c(~σk) as the local upscaling criterion. Inthis case, our formulation is a type of extended upscaling procedure within theclassification of upscaling techniques proposed in [49, 56].Once again, the parameter estimation problem (3.17) satisfies that the numberof parameters is less than (or equal) to the number of data and it is a well-posedleast-squares problem, hence no additional regularization is required.In order to solve the discrete optimization problem (3.17) using a gradientbased optimization method (see [107, 117, 137] for details), we first compute thegradient and an approximation to the Hessian of the local upscaling criterion c asdiscussed in [68]. The gradient of the local upscaling criterion (3.17) is given by∇c(~σk) =12∑l=1Real(J∗l (~σk)[dl(ΣHk (~σk))−dl(Σhk)]) , (3.18)where Jl represents the lth sensitivity matrix defined as the Jacobian resultingafter deriving the lth residual, dl(ΣHk (~σk))− dl(Σhk), with respect to ~σk, and ∗denotes the conjugate transpose operator. Since each sensitivity matrix has di-mensions 12 × number of data (e.g. 12 × 6 or 12 × 12), we can compute themexplicitly.The Gauss-Newton approximation to the Hessian of the local upscaling crite-rion (3.17) is given byHc(~σk)≈12∑l=1Real(J∗l (~σk)Jl(~σk)). (3.19)where Jl is the lth sensitivity matrix defined as before.The pseudo-code in Algorithm 1 summarizes the steps to compute an up-scaled conductivity ΣHk in the coarse cell ΩHk when a discrete fine-mesh conduc-tivity modelΣh is given. Observe that all the calculations are done inΩH,extk , whereeach optimization problem is small and can be solved quickly. Furthermore, sincethe problem defined for each coarse cell is independent, the upscaling procedure62Algorithm 1 Computation of the upscaled conductivity ΣHk in the coarse cell ΩHk :1: Choose the size of the extended local domain ΩH,extk (Figure 3.4) where thecoarse cell, ΩHk , is embedded.2: Compute sets of discrete fine and coarse EM responses using the fine-meshconductivity in ΩH,extk . To do so, forward model the twelve local Maxwell’sproblems (3.13)-(3.15) defined onΩH,extk using the MFV method as outlined inSection 2.4. This yields the sets of discrete electric fields (Ek = {ek1 , . . . ,ek12})and magnetic fluxes (Bk = {bk1 , . . . ,bk12}).3: Choose the type of data to be matched in the local upscaling criterion (3.17)according to the context of the given simulation. That is, chose either the setof total electric fields as defined in (3.16a), or the set of total magnetic fluxesas defined in (3.16b). Compute the fine and coarse data to be matched usingthe discrete fields and fluxes obtained in the previous step.4: Optimize the constrained parameter estimation problem (3.17) to obtain thedesired anisotropic upscaled conductivity ΣHk . Depending on the choice of theoptimization method to solve this problem, this step may require the computa-tion of the gradient (3.18), and/or the Hessian approximation (3.19). Note thatsolving such an optimization problem involves performing steps 2 and 3 (forthe upscaled conductivity model ΣHk ) to compute dl(ΣHk)at each iteration.can be done in parallel.Once all the upscaled conductivities are computed, we assemble the coarse-mesh conductivity model. Using this coarse conductivity model, the original prob-lem can then be discretized and solved on the coarse mesh using a traditionaldiscretization approach, such as the MFV method (Section 2.4).Remarks1. Observe that the upscaled conductivity constructed using formulation (3.17)can result in a SPD diagonal matrix or a full matrix depending on the prob-lem. That is, the formulation leads to the appropriate form of the best up-scaled conductivity (in the sense of equation (3.17)) for the problem at hand.This feature of the proposed upscaling formulation can be advantageous asfor some averaged-based upscaling methods proposed for flow in porousmedia problems this is not the case. For example, two-point flux approxi-mation upscaling methods assume that the upscaled quantity is a diagonal63matrix. This assumption may not be sufficient to obtain an accurate so-lution in the situations where a full-tensor is more appropriate. For thosecases, full-tensor effect or full anisotropy of the upscaled permeability canbe achieved by post-processing the diagonal tensor using fine-mesh infor-mation where care must be taken to preserve the symmetry of the matrix[35, 36, 49, 56].2. The proposed upscaling formulation (3.17) can be extended to constructcomplex anisotropic conductivities as well. One way to do so is by adaptingthe parametrization of the upscaled conductivity given by (3.12) accordinglyand solving the resulting optimization formulation in a similar manner asdescribed before. That is, rather than using six scalars to parametrize ΣHk ,we use twelve scalars where six of them are used to define the real partof ΣHk and the other six scalars define the imaginary part of ΣHk . Usingsuch parametrization leads to a real optimization problem with a similar op-timization formulation as before where ~σk ∈ R12 and Real(ΣHk (~σk)) shouldbe SPD. The computation of the discrete data to be matched in the upscal-ing criterion can be done in a similar manner as described before; however,the interested reader can find more detailed information on how to computeEM fields when the electrical conductivity is a complex quantity in [68].3. The data to be matched by the local upscaling criterion (3.17) can also bechosen among the total magnetic fields, or the total electric fluxes on ∂ΩHk .The total magnetic field can be computed as the integral of the magneticfield, ~H, over the twelve edges of ∂ΩHk (i.e., replace ~E with ~H in (3.16a)).The total electric flux can be computed as the integral of the electrical cur-rent density, ~J, over the six faces of ∂ΩHk (i.e., replace ~B with ~J in (3.16b)).The discretization of the magnetic field and the electrical current density canbe done using the ~H-~J formulation of the quasi-static Maxwell’s equationsand the MFV in a similar way as described in Section 2.4 for the ~E-~B for-mulation. In Section 3.5.1, we show two examples that use total magneticand electric fluxes to compute upscaled electrical conductivities on a singlecoarse cell.643.5 Numerical Results in 3DIn this section, we demonstrate the performance of the local upscaling frameworkproposed in Section 3.4 by simulating EM responses for two 3D settings. Thefirst setting illustrates the use of the framework to upscale fine-scale conductiv-ity structures on a single coarse-mesh cell, and the impact that the choice of theupscaling criterion has on the resulting upscaled quantity. The second settingdemonstrates how the framework can be combined with an adaptive mesh refine-ment technique to construct anisotropic, coarse-mesh conductivity models for alarge-loop EM survey of a mineral deposit on a geologically-rich background.3.5.1 Simulations on a Single Coarse-mesh CellWe demonstrate the upscaling procedure and the impact of the choice in upscal-ing criterion on a single coarse-mesh cell. To do so, we use two electrical conduc-tivity models: one of an isolated conductive block on a resistive background, andone of a conductive sheet on a resistive background. These conductivity modelsare visualized in Figure 3.6 (a) and (b), respectively.In both examples, the local upscaling domain is given by a cuboid coarse-mesh cell with dimensions (100 m)3. We denote the local upscaling domain asΩHk . We extend the local domain a further 50 m along each direction such thatΩHkis positioned in the center (Figure 3.6). This conforms the local extended domain,which we denote as ΩH,extk . Note that the local extended domain has dimensions(200 m)3.We discretize the local extended domain using uniform fine-mesh cells, whichare (12.5 m)3. Hence, the total number of fine cells in ΩH,extk is 163, and in ΩHkis 83. As discussed in Section 3.4, we assume the fine-scale conductivity modelto be isotropic and the upscaled conductivity to be anisotropic in the coarse cellwe aim to upscale. The magnetic permeability takes its free space value, namelyµ = 4pi×10−7 Vs/Am.To construct an anisotropic upscaled conductivity using the upscaling criterion(3.17), we must first select the data to be matched on it. For the following twoexamples, we consider two types of data: the total magnetic flux on each of thesix faces of ΩHk , and the total current density flux on each of the six faces of ΩHk .65Figure 3.6: Cross sections of the fine scale conductivity model for (a) theisolated block and (b) the sheet. Conductive bodies are displayed inred. The resistive background is displayed in blue. The coarse-meshcell, ΩHk , which we aim to upscale is outlined in white.The total flux (either current density or magnetic field) over each face is definedas the surface integral of the flux through that face (see (3.16b)). Throughout thissection, we say that we use the b−criterion when the data to be matched in theupscaling criterion (3.17) are total magnetic flux. Similarly, we say that we use thej−criterion when the data to be matched in the upscaling criterion (3.17) are totalcurrent density flux.Once we selected the necessary elements to setup the local upscaling frame-work proposed in Section 3.4, let us see how the framework performs for the twoconductivity models we are interested in.Conductive Block in a Resistive BackgroundFirst, we study the case of an isolated conductive block in a resistive background.To setup this example, we assume that at the center of ΩHk there is a (50 m)3conductive block, as shown in Figure 3.6(a). The surrounding material is resistive(10−4 S/m). We use a frequency of 1 Hz and consider three different electricalconductivities of the block: 10−2, 10−1 and 1 S/m.Following Algorithm 1, next we need to solve the optimization problem (3.17)66to construct anisotropic upscaled conductivities, as defined in equation (3.12), forthe b−criterion and j−criterion, and each of the three conductivity values of theblock, respectively. The initial conductivity model, ~σ0, is given by a 3×3 diagonalmatrix with all its main diagonal entries being equal to the geometric mean of thefine-mesh conductivity in ΩHk . We choose the geometric mean as it is often usedto approximate the homogenized permeability for fluid flow problems with randomheterogeneous media with reasonable accuracy [163].In order to solve the constrained optimization problem (3.17), we tested botha projected steepest descent and a projected Gauss-Newton method (see [68,137]). Both methods were implemented with a backtracking line search. Theprojection in this case consisted of computing the eigenvalues of the 3×3 matrixdefined by (3.12) after ~σk was computed at each iteration to ensure the positivedefiniteness of such matrix. When necessary, the nonpositive eigenvalues of suchmatrix were modified to be positive. For both optimization methods, the iterationwas halted when‖~σk−~σk+1‖2 < 10−8 or ‖∇c(~σk)‖2 < 10−8‖∇c(~σ0)‖2. (3.20)Although the use of the projected Gauss-Newton method lacks theory for conver-gence in this case, in practice we observe that this method reaches results similarto those obtained with the projected steepest descent method but in a much fasterway. For this reason, we decided to use the projected Gauss-Newton method tocarry out the simulations.For each of the three conductivity values considered, the resulting upscaledconductivities are diagonal matrices with all diagonal entries being equal. Thisresult is expected due to the symmetry of the fine-mesh conductivity model. Whenusing the b−criterion, the stopping criterion (3.20) was fulfilled for each of thethree conductivity values of the block considered after 18, 4 and 26 iterations,respectively. When using the j−criterion, the stopping criterion (3.20) was fulfilledfor each of the three conductivity values of the block considered after 17, 14 and20 iterations, respectively. The values of the main diagonal entries of each ofupscaled conductivity are shown in Figure 3.7.Clearly, there is a large discrepancy between the upscaled conductivities con-67Figure 3.7: Upscaled conductivity results found using the b and j criteria forthe conductive block model. Note that the upscaled conductivity in thiscase can be defined by a scalar, as the matrix recovered was diagonal,with all diagonal elements being equal.structed using the b− and j− upscaling criteria for each of the three block con-ductivities. Physically, this discrepancy can be reconciled by recognizing that thecurrent density and magnetic flux density reflect different physical processes. Cur-rent is the flow of electrical charges through a material. In the case of a conductiveblock in a resistive background, the current must be driven through the resistivebackground irrespective of the conductivity of the block. Thus, the resistive back-ground dominates the upscaled conductivity constructed using the j-criterion. Onthe other hand, magnetic flux is produced as a result of induced currents. Cur-rents can be induced in the conductive block regardless of the conductivity of thebackground. Therefore, we see that the conductivity of the block dominates theupscaled conductivity recovered using the b-criterion.Conductive Sheet in a Resistive BackgroundNext, we examine the case of a conductive sheet in a resistive background posi-tioned at the center of ΩHk , as shown in Figure 3.6(b). The sheet has dimensions150 m × 150 m × 50 m. It is stopped short of the boundary of the coarse cell68in order to avoid applying the boundary conditions directly on the sheet. For thissimulation, we also use a frequency of 1 Hz.Now that the model is no longer identical in each direction, we expect to con-struct an anisotropic upscaled conductivity, as defined in (3.12), where σ1 = σ2which is distinct from σ3.Typically, one would draw the analogy between layered conductors and a sim-ple circuit model. In this case, we would expect that the components of the up-scaled conductivity matrix, σ1 and σ2 would conform closely the approximationof conductors in parallel, while σ3 would be similar to the approximation of con-ductors in series. Note however, that this approximation accounts for only onephysical behavior: galvanic current flow, to which the j-criterion is sensitive. Itdoes not account for any inductive currents, to which the b-criterion is sensitive.To investigate, we again assign the electrical conductivity of the backgroundto be 10−4 S/m and examine three conductivities for the sheet: 10−2, 10−1 and1 S/m. Using both the b− and j− criteria, we perform the upscaling process asindicated in Algorithm 1, and recover the upscaled conductivities shown in Figure3.8.To obtain the results described above we use the stopping criterion (3.20) andthe same initial conductivity model (~σ0) as before. We also run the simulations us-ing both projected steepest descent and projected Gauss-Newton methods. Onceagain, we notice that the results obtained with projected Gauss-Newton are simi-lar to those obtained with projected steepest descent, but the convergence whenusing projected Gauss-Newton is faster. When using the b−criterion, the stop-ping criterion (3.20) was fulfilled for each of the three conductivity values of thesheet considered after 24, 54 and 142 iterations, respectively. When using thej−criterion, the stopping criterion (3.20) was fulfilled for each of the three conduc-tivity values of the sheet considered after 29, 67 and 181 iterations, respectively.As expected, the resulting upscaled conductivity is a diagonal SPD matrixwith two unique entries: σ1 = σ2 and σ3 (see (3.12)). Figure 3.8(a) shows theσ1 = σ2 entries of the matrix, and Figure 3.8(b) shows the σ3 entry. The valuesconstructed using the j-criterion conform well to the parallel and series circuitapproximations for σ1 = σ2 and σ3, respectively. For each scenario the value forσ3 constructed using the j-criterion is nearly identical to the resistive background,69(a)(b)Figure 3.8: Upscaled conductivity results found using the j and b criteria forthe conductive sheet. Since the upscaled conductivity is a diagonal3×3 SPD matrix, it can be described by two positive scalars: σ1 = σ2and σ3. σ1 = σ2 is shown in plot (a), and σ3 is shown in plot (b).as is to be expected using a series circuit approximation. Similar to the blockmodel results, this is as a consequence of having to drive the current through theresistive background in the z direction, regardless of the conductivity of the sheet.On the other hand, for the entries σ1 and σ2, the conductivity of the sheet has a70large impact on the value we construct. In this case, the conductive sheet forms aconnected pathway from one side of the cell we aim to upscale to the other alongboth horizontal directions. As a result, current is channeled along this pathway,causing the conductivity of the sheet to have a large impact on the σ1 and σ2entries of the upscaled conductivity, as shown in Figure 3.8.Using the b-criterion, we see that the conductivity of the sheet has a significantimpact on both the σ1=σ2 (Figure 3.8(a)) entries and the σ3 entry (Figure 3.8(b)),contrary to the parallel-series circuit approximations. This is again because themagnetic flux density is sensitive to inductive currents. Since the sheet has a finitethickness, these currents can be induced in any direction, and therefore contributeto the upscaled values σ1, σ2 and σ3 we construct using the b-criterion.Lessons LearnedAs discussed in the block and sheet examples, the current density and magneticflux density each sample the conductivity structure differently, and are thereforesensitive to different features of the fine-mesh conductivity model. As a result, theconstructed upscaled conductivities using either upscaling criteria may vary overorders of magnitude. Therefore, for a given fine-scale conductivity structure, thereis not a unique upscaled conductivity which completely describes it.3.5.2 Simulations Using the Synthetic Lalor ModelIn this section, we demonstrate how to use the upscaling framework introducedin Section 3.4 in combination with an adaptive mesh refinement technique to con-struct anisotropic, coarse-mesh electrical conductivity models from a given fine-scale conductivity model of a mineral deposit. We use these coarse-mesh con-ductivity models to predict magnetic field data for a large-loop EM survey. Sinceno analytical solutions are available for this example, we compare our results withthose obtained on a fine mesh. We also compare our results with those obtainedfrom applying simple average-based upscaling approaches.We construct a synthetic electrical conductivity model based on the inversionresults of EM field measurements over the Canadian Lalor mine obtained by [177].The Lalor mine targets a large zinc-gold-copper deposit that has been the subject71Figure 3.9: Subsurface part of the synthetic electrical conductivity Lalormodel and large-loop EM survey setup. The model is discretized ona fine OcTree mesh with 546,295 cells. The conductivity varies overfive orders of magnitude throughout the whole model.of several EM surveys. The synthetic conductivity model, shown in Figure 3.9,has an area with non-flat topography and extends from 0 to 6.5 km along the x, yand z directions, respectively. The model comprises air and the subsurface thatis composed of 35 geologic units. The unit with the largest conductivity valuerepresents the mineral deposit, which is composed of three bodies. We assumea conductivity of 10−8 S/m in the air. The subsurface conductivity values rangefrom 1.96× 10−5 to 0.28 S/m. We assume that the magnetic permeability takesits free space value, namely µ = 4pi×10−7 Vs/Am.We consider a frequency-domain, large-loop EM survey, where we use a rect-angular transmitter loop with dimensions 2 km × 3 km, operating at the frequen-cies of 1 and 20 Hz. The transmitter is placed on the earth’s surface and it is72centered above the largest body of the mineral deposit, as shown in Figure 3.9.Inside the loop, we place a uniform grid of receivers that measure the three com-ponents of the magnetic flux (~B = [Bx,By,Bz]>). The receivers are separated by50 m along the x and y directions, respectively. To reduce the effect of the im-posed natural boundary conditions (2.5), we embed the survey area into a muchlarger computational domain, which replaces the true decay of the fields towardsinfinity (Figure 3.9).We discretize the synthetic Lalor electrical conductivity model using a stag-gered fine OcTree mesh, which allows for an adaptive local refinement of the meshwhere the electrical conductivity and the EM responses vary drastically. We basethe estimate for the proper cell sizes of the mesh on the skin depth [169]. Practi-cal experience on mesh design for EM problems reported in [68] suggests that thesmallest cell size in the mesh should be a quarter the minimum skin depth. Weconsider the largest background conductivity value (4.5×10−3 S/m) and calculateskin depths of 7,498 and 1,677 m for the frequencies of 1 and 20 Hz, respectively.Hence, we use cells of size (50 m)3 within the survey area and at the interfacesof the model where the conductivity varies, the rest of the domain is padded withgradually increasing OcTree cells. The fine OcTree mesh is illustrated in Figure3.9. This mesh has 546,295 cells.We aim to estimate the secondary magnetic flux induced by the ground in thesurvey area. For this purpose, we simulate two sets of the magnetic flux datafor each frequency. The first data set considers the conductivity model includingall geologic units, and the second data set considers the conductivity model of ahalfspace with an earth conductivity of 10−2 S/m. Each of these two data setsconsists of the measurements of ~B taken at the receiver locations. The secondarymagnetic flux induced by the ground at the survey area, denoted as ∆~B, is thencomputed by subtracting the two data sets.In order to construct a single anisotropic, coarse-mesh conductivity modelusing the upscaling methodology introduced in Section 3.4, we need to choosethe following parameters: (a) a suitable coarse mesh, (b) the size of the extendedlocal domain, and (c) the data to be matched in the upscaling criterion (3.17).As a coarse mesh, we consider an OcTree mesh that is nested within thefine OcTree mesh previously described. Since we are interested in accurately73Figure 3.10: Subsurface part of the synthetic electrical conductivity Lalormodel and large-loop EM survey setup. The model is discretized on acoarse OcTree mesh with 60,656 cells. The conductivity varies overeight orders of magnitude throughout the whole model.simulating magnetic flux data in the survey area, the OcTree mesh is designed tomaintain the fine mesh resolution of (50 m)3 inside the area of interest, whereasthe rest of the domain is filled with increasingly coarser cells. In total it containsonly 60,656 cells — this number of cells is roughly 10% of the number of cellsin the fine OcTree mesh. Figure 3.10 shows the coarse OcTree mesh. Note thatthe coarse mesh is not refined outside the survey area where a large conductivitycontrast is present in the model. For example, at the interface between the highlyconductive gold units and the more resistive background, and at the air-earthinterface which contains a not flat area. These interfaces are not represented inthe coarse OcTree mesh. We challenge the upscaling procedure with the largeconductivity contrast across the subsurface interfaces.74Next, we need to choose the size of the extended domain to solve the localMaxwell problems, that is, the number of fine-mesh padding cells by which we ex-tend every coarse cell to be upscaled as shown in Figure 3.4. To investigate theeffect of this size on the resulting upscaled conductivity, we performed the exten-sion using the sizes of 4 and 8 padding cells. These local extended domain sizesextend each (200 m)3-coarse cell by one and two coarse cells, respectively. The(200 m)3 cells are the majority of the coarse cells where the largest conductivitycontrast happens in this setting (Figure 3.10).Since the receivers measure magnetic flux data, we consider the total mag-netic flux as defined in (3.16b) through each of the faces of the coarse cell tobe upscaled as the data to be matched by the upscaling criterion (3.17). By do-ing so, we have connected the upscaling criterion to the large-loop EM surveyconfiguration.Once the necessary elements to setup the upscaling framework are selected,the construction of the desired coarse-mesh conductivity models is completedby solving the constrained parameter estimation problem (3.17) for each of thecoarse cells separately, each individual size of the extended domain, and eachindividual frequency. Each individual parameter estimation problem is solved asindicated in Algorithm 1. To solve each local parameter estimation problem weuse the same stopping criterion and optimization parameters as described in theprevious section.Next, we use the coarse conductivity models obtained previously to estimate∆~B. To do so, we apply the MFV discretization method on the described coarseOcTree mesh as outlined in [72, 86]. This discretization yields linear systems ofequations with 169,892 edge degrees of freedom (DOF), which we solve usingthe parallel sparse direct solver MUMPS [3]. The total average run times persingle simulation for extended domain sizes of 4 and 8 padding cells are 5.3 hand 3.7 days, respectively, on a two hexa-core Intel Xeon X5660 CPUs at 2.8 Hzwith 64 GB shared RAM using the parallel computing toolbox of MATLAB ([162]).The right-hand panel of Figure 3.11 shows the magnitude of ∆~B obtained for thefrequencies of 1 and 20 Hz for the extended domain size of 8 padding cells.To evaluate the accuracy of our results, we compute a reference (fine-mesh)solution. To do so, we apply the MFV method on the described fine OcTree mesh75(Figure 3.9) as discussed in [72, 86]. This yields linear systems of equationswith roughly 1.5 millions DOF, which we also solve using MUMPS. The averagecomputation time per single simulation is 704.5 s on the same machine. The left-hand panel of Figure 3.11 shows the magnitude of ∆~B obtained for the frequenciesof 1 and 20 Hz.We also carry out MFV forward modeling simulations using homogenizedelectrical conductivity models that we construct using volume-arithmetic, -geometricand -harmonic averaging of the fine-mesh conductivity inside each coarse cell ofthe OcTree mesh shown in Figure 3.10. The total average run time per singlesimulation is 125.5 s on the same machine.Table 3.4 shows the relative errors in Euclidean norm for the magnitude of∆~B obtained from comparing the reference (fine-mesh) solution with the differenthomogenized solutions for each frequency and local extended domain size. Therelative error is computed as the ratio of the norm of the difference in magnitudeof the fine- and coarse-mesh data to the norm of the fine-mesh data in magnitude.Table 3.4: Relative errors in Euclidean norm for the magnitude of the sec-ondary magnetic fluxes induced by the ground in the survey area, ∆~B.The relative errors are given in per cent.Conductivitymodel1 Hz 20 HzArithmetic 9.028 1.534Geometric 8.990 0.594Harmonic 8.987 0.403Upscaling (4 padding cells) 8.989 0.532Upscaling (8 padding cells) 8.991 0.383From Table 3.4 we see that while upscaling has an effect when using 20 Hzdata, the effect is smaller when considering the 1 Hz data. This should not comeas a surprise, as the fields at 1 Hz are mainly in the real component, and areless sensitive to fine-scale variations in conductivity than at 20 Hz. Indeed, for theMagnetostatic case, the magnetic fields are conductivity independent [68]. As aresult, the error that we observe in this frequency is mostly due to discretizationerror. However, when considering 20 Hz, the effect of using an appropriate aver-76aging scheme is more evident and in fact, our averaging scheme does better thanother averaging schemes. Nonetheless, for this case, it is surprising to see howwell simple harmonic averaging did. We suspect that, unlike the 1D case whereharmonic averaging performs poorly, the 3D conductivity model requires averag-ing over a much smaller, local area which leads to a much smaller difference inthe data.77Reference solution Upscaled solution (8 padding cells)1Hz20HzFigure 3.11: Magnitude of the secondary magnetic flux induced by the ground in the survey area, ∆, for the large-loopEM survey configuration. First and second rows show the results for 1 and 20 Hz, respectively. The left-handpanel shows the reference solution computed using the MFV method on the fine OcTree mesh with 546,295cells. The right-hand panel shows the upscaled solution computed using 8 padding cells on the coarse OcTreemesh with 60,656 cells.783.6 DiscussionThe global and local upscaling formulations proposed in this chapter are appro-priate in different situations. For example, Section 3.3 shows the benefit of us-ing the global upscaling formulation (in 1D) when upscaling well log conductivitydata. On the other hand, Section 3.5.1 shows the benefit of using the local up-scaling formulation to construct 3D anisotropic coarse-mesh conductivity modelsfor a large-loop EM survey of a mineral deposit. This investigation supports theremarks made in the field of fluid flow problems regarding the use of upscalingmethods (cf. [49]): the upscaling formulation to use on a particular problem de-pends on the simulation question being addressed and the level of detail that canbe accommodated in the coarse conductivity model.The 1D and 3D experiments presented show that the coarse-mesh modelsconstructed with the proposed upscaling framework yield accurate approximationsto the EM responses that are comparable to those obtained by using MFV on afine mesh in the forward modeling process, and that the size of the problem can bereduced significantly, specially when upscaling is combined with an adaptive meshrefinement technique, such as OcTree. For the examples presented, the size ofthe coarse-mesh system solved was roughly 10% of the fine-mesh system size,while the relative errors (in the secondary fields) were less than 5%. That is, thecoarse conductivity models are able to emulate the behavior of the heterogeneitypresent in the prescribed fine-mesh electrical conductivity model.The 3D local upscaling framework has two main issues. First, it constructs anupscaled conductivity that depends significantly on the set of boundary conditionsimposed to compute the synthetic data used in the upscaling criterion. This issueis consistent with what has been reported in the upscaling procedures developedin the community that simulates flow in porous media [49, 56, 175]. I use the set ofstandard bilinear decaying functions on a coarse cell as the set of boundary con-ditions based on the arguments described in Section 3.4. I recognize that such aset of boundary conditions may not be the most appropriate for constructing accu-rate coarse models for all cases. However, for the experiments presented, theseboundary conditions give reasonable estimates. Second, the proposed upscalingmethod is more computationally expensive than simple average-based upscaling79methods, as one solves a local optimization problem in each coarse cell. How-ever, since each local problem is formulated independently of the others, one canreduce the cost by using a more efficient parallel implementation of the methodon a more powerful machine.3.7 SummaryThis chapter proposes a least-squares formulation for the upscaling problem anddevelops a numerical framework to construct accurate coarse-mesh electricalconductivity models based on prescribed fine-mesh ones for a broad range ofquasi-static EM geophysical problems in the frequency domain. In practice, sim-ulating these types of problems is computationally expensive; they often considerhighly heterogeneous geologic media that require a very large and fine mesh tobe discretized accurately.In the proposed framework, we pose upscaling as a parameter estimationproblem. Thus, a coarse-mesh electrical conductivity model is obtained by solv-ing an optimization problem for each coarse-mesh cell. The optimization crite-rion (i.e., upscaling criterion) can be customized to construct isotropic or fullyanisotropic real or even complex upscaled quantities to approximate any of theEM fields and/or fluxes depending on the geophysical EM experiment of inter-est. The computation of the upscaled conductivity can be done in a global or anextended setting. Contrary to other least-squares upscaling formulations in the lit-erature for fluid flow problems, the formulation proposed in this investigation doesnot require regularization and it is practical to tackle 3D geophysical EM problemswith arbitrary electrical conductivity structures.The upscaling framework demonstrates that the construction of upscaled quan-tities should be specific to, and highly dependent on, the purpose of the simula-tion and the EM experiment configuration. In fact, this investigation shows thatdifferent EM experiments use different upscaling criteria that result in different up-scaled quantities. The proposed upscaling formulation is a reflection of the factthat an upscaled quantity is a property that one constructs, and it will take on dif-ferent values, depending on how one formulates the problem. The choices of theEM responses of interest, boundary conditions, and type of the upscaled quan-80tity employed (i.e., isotropic or anisotropic) all influence the nature of the resultingupscaled conductivity model. As a result, for a given fine-scale conductivity struc-ture, there is no unique upscaled model which completely describes it. I believethat the proposed framework can be used to tackle upscaling problems with adifferent perspective when care is taken to properly define the upscaled quantityneeded.This investigation demonstrates that upscaling methods can be used for solv-ing frequency-dependent, quasi-static EM geophysical problems with highly dis-continuous electrical conductivity, especially given the work done on improvingthe computation performance and accuracy of the proposed framework. Simul-taneously, Haber and Ruthotto ([73]) showed that multiscale techniques are alsoa feasible method to tackle this problem with some additional advantages. Multi-scale methods are further investigated in the next chapter.81Chapter 4A Multiscale Finite VolumeMethod with OversamplingEvery pawn is a potential queen. — James Mason4.1 OverviewThis chapter1 proposes a multiscale FV method with oversampling for faster sim-ulation of the quasi-static EM responses on a coarse mesh.The main contribution of this chapter is to show how the core mathemati-cal concepts used to develop oversampling techniques for flow in porous mediaproblems can be extended to geophysical EM problems, as well as to show howmultiscale techniques can be combined with OcTree in a parallel environment totackle more challenging geophysical EM simulations.This chapter starts by providing a summary of the multiscale FV method pro-posed by Haber and Ruthotto [73] for geophysical EM problems. Then, the pro-posed oversampling technique for such method is developed. Finally, the perfor-mance of the multiscale FV with oversampling method is demonstrated by usingtwo 3D synthetic electrical conductivity models: a) the Lalor model, and b) one1This chapter contains extended and revised versions of the material published in [30], [31] and[33]. I am the lead author in all of these publications.82with a random isotropic heterogeneous medium. Both examples show the fea-sibility of combining the proposed multiscale FV with oversampling method withOcTree meshes in a parallel environment to boost its performance.4.2 An Overview of the Multiscale Finite Volume MethodHaber and Ruthotto ([73]) adapted the general lines proposed by Hou and Wu([89]), Jenny et al., ([101]), and MacLachlan and Moulton ([123]), where multiscaleFE and FV methods are developed for elliptic problems with strongly discontinu-ous coefficients, to develop a multiscale finite volume method (MSFV) that fits thestaggered discretization of vector fields typically used in the MFV discretizationmethod introduced in Section 2.4. Since we use the MSFV method as a buildingblock for the oversampling technique we propose in this work, we provide nextan overview of this method, which can be summarized in the following four steps(Figure 4.1).Figure 4.1: Schematic representation of the procedure to implement theMSFV method.83First, let us assume a coarse mesh, MH , nested into a fine mesh, Mh, i.e.,MH ⊂Mh (Step 1 in Figure 4.1). Mh accurately discretizes the features in themodel where the electrical conductivity varies, and MH is a user-chosen meshthat is typically much coarser than the fine mesh and satisfies the guidelines formesh design in the areas where the EM responses are measured. In particular,MH =∪Nk=1ΩHk , where N is the number of coarse-mesh cells and ΩHk denotes thekth coarse-mesh cell;Mh =∪ni=1Ωhi , where n is the number of fine-mesh cells andΩhi denotes the ith fine-mesh cell; and N n. The MSFV method was originallydeveloped for nested tensor meshes, we will show an example in Section 4.4where we use nested OcTree meshes as the mesh setup.Second, for each coarse-mesh cell, ΩHk ; k = 1, . . . ,N, we independently solvea local version of the source-free Maxwell system subject to a set of twelve lin-early independent, non-homogeneous Dirichlet boundary conditions (one for ev-ery edge of ΩHk ) given by∇× ~Ekl + ıω~Bkl = ~0, in ΩHk , (4.1)∇× (µ−1,k~Bkl )−Σk~Ekl = ~0, in ΩHk , (4.2)~Ekl (~x)×~n = ~Φl(~x)×~n, ∀~x ∈ ∂ΩHk , (4.3)with l = 1, . . . ,12. Here, ~x = (x1,x2,x3), ∂ΩHk denotes the boundary of ΩHk ,~Bkl , ~Ekl , µk, Σk, ω , ı and~n are defined as before in Section 2.3; and ~Φl is a vectorfunction that denotes the lth Dirichlet boundary condition as defined in Table 3.3.This set of boundary conditions form the natural basis functions for linear edgedegrees of freedom of a hexahedral finite element [104, 131]; therefore, they canbe used to model general normal-linearly varying EM responses. From Figure3.5, we observe that each ~Φl takes the value 1 along the tangential direction tothe lth edge of ΩHk and decays linearly to 0 in the normal directions to the sameedge.To numerically solve the twelve local Maxwell systems (4.1)-(4.3), we forwardmodel them using the fine mesh contained in ΩHk and the MFV method as dis-cussed in Section 2.4 (Step 2 in Figure 4.1). That is, the set of discrete solutionsfor the electric field, {ek1, . . . ,ek12}, can be obtained by solving twelve linear sys-84tems of the form (2.27), where each ekl ; l = 1, ...,12, is a vector whose lengthequals the number of fine-mesh edges in ΩHk . We use the MFV method be-cause it provides a consistent and stable discretization of the Maxwell’s equationswith highly discontinuous coefficients. However, other traditional edge-based dis-cretization methods, such as FE ([104, 131]), can be used for this part as wellas long as they provide a consistent, conservative and stable discretization forour Maxwell system. In this chapter, we refer to the set of discrete solutions{ek1, . . . ,ek12} as the multiscale basis functions for the cell ΩHk . The MSFV methodproposed by Haber and Ruthotto ([73]) only covers the necessary steps to com-pute multiscale basis functions for the electric field. They do not compute basisfunctions for the magnetic flux, and neither do we. To see full derivation details ofthese multiscale basis functions, we refer the interested reader to the body of thepaper [73].Note that using formulation (4.1)-(4.3) implies that ~Ekl is oscillatory at the in-terior of the coarse cell ΩHk , and that it coincides with the natural basis functions{~Φl}12l=1 at the boundary of ΩHk , that is~Ekl ·~τedgem = δlm; l,m = 1, . . . ,12, (4.4)where ~τedgem is the unit tangent vector to the mth edge of ∂ΩHk , and δlm is theKronecker delta that takes the value 1 when l = m and 0 otherwise. Naturally, themultiscale basis functions also satisfy these properties. It follows that the tangen-tial components of the multiscale basis functions are continuous at the boundariesof the coarse-mesh cells.As shown in [51, 69, 89], the multiscale basis functions can be arranged asthe columns of a local coarse-to-fine interpolation matrix Pk, that isPk =[ek1, . . . ,ek12], (4.5)for the fine-mesh electric field in ΩHk (Step 2 in Figure 4.1). This type of interpo-lation is also known as operator-induced interpolation, which was originally devel-oped in the Multigrid community for the diffusion equation with strongly discontin-uous coefficients [2, 44].85Once we have computed a local interpolation matrix for each coarse cell (4.5),the third step is to assemble a global coarse-to-fine interpolation matrix, P (Step 3in Figure 4.1). The continuity of the tangential components of the multiscale basisfunctions at the boundaries of the coarse cells is a necessary requirement for theproper assembly of P in this step [51].We point out that the calculations involved to compute the local interpolationmatrices (Pk; k= 1, . . . ,N) as defined in (4.5), are done locally inside each coarsecell independently of each other, hence they can perfectly be done in parallel. Thisgreatly reduces the overhead time in constructing each Pk in practice.The fourth step is to use the global interpolation matrix P as a projection matrixwithin a Galerkin approach ([51]) to construct a coarse-mesh version of the fine-mesh system (2.25) that is much cheaper to solve as followsAHeH =(P>Ah(Σh)P)eH = P>qh. (4.6)The superscripts H and h denote dependency to the coarse and fine meshes,respectively, the vector qh and the system matrix Ah(Σh) are defined as in (2.25),and eH denotes the coarse-mesh electric field. To compute the coarse-meshmagnetic flux, bH , we use eH in (2.23).As shown in [69], the fine-mesh electric field, eh, can be obtained from thesolution to the coarse-mesh problem as followseh = PeH . (4.7)To compute the fine-mesh magnetic flux, bh, we use eh in (2.23). This concludesthe overview of the MSFV method proposed by Haber and Ruthotto in [73].Note that we use the exact same formulation for the local Maxwell’s problems(4.1)-(4.3) for both the MSFV method and the local upscaling method in 3D in-troduced in Section 3.4 (see equations (3.13)-(3.15)). Although the formulationof these local problems is the same, they differ in the context in which they wereformulated. The difference is that in the local upscaling method the formulationcomes from using physical arguments that describe the behavior of our problem todefine the elements needed in the upscaling criterion (3.17); whereas in the MSFVmethod the formulation comes from the Galerkin framework we use to compute a86set of basis functions that are necessary to construct the projected system (4.6).In the above process, we opt for the construction of the operator-induced in-terpolation matrix P; however, it is possible to avoid its construction. As shownin Chapter 2 of the book by Efendiev and Hou ([51]), it is possible to assembledirectly the matrix AH of the coarse-mesh system in (4.6) through locally pro-jected stiffness matrices, which are generated using Pk. The former approach ispreferred when the solution is needed in the coarse mesh only as it reduces thestorage requirements. In Section 4.4, I show examples of an application of thistype. If the fine-mesh solution is needed, then the first approach is preferred.The accuracy of the solutions obtained with multiscale FV/FE methods de-pends on the choice of boundary conditions used to construct the multiscale ba-sis functions (second step outlined before) for each coarse cell. If these boundaryconditions fail to reflect the effect of the underlying media heterogeneity containedby the coarse cell on the physical responses, multiscale procedures can havelarge errors [51, 89].Researchers in the field of multiscale methods for elliptic problems have notedthat by choosing a set of linear boundary conditions for the construction of the mul-tiscale basis functions, a mismatch between the exact solution and the discretesolution across the coarse-cell boundary may be created, thus yielding to inac-curate solutions. The error analyses presented in [89] and in [51] demonstratethat the source of inaccuracy in the solution comes from resonance errors; thatis, errors that appear when the mesh size and the wavelength of the small-scaleoscillation of the media heterogeneity are similar. A solution in such cases is touse oversampling techniques for the construction of the multiscale basis functions[51, 79, 88, 89, 101].In the next section, we discuss the case where the choice of linear boundaryconditions to construct the multiscale basis functions may yield inaccurate solu-tions using the MSFV method, and we develop an oversampling technique to fixthis accuracy issue.874.3 The Oversampling MethodAs discussed in the previous section, the MSFV method imposes linear bound-ary conditions to the local Maxwell’s formulation (4.1)-(4.3) used to compute themultiscale basis functions inside each coarse-mesh cell. Note that by choosinglinear boundary conditions for the multiscale basis functions, the MSFV methodassumes that the tangential components of the electric field behave linearly atthe interfaces between coarse-mesh cells. However, this assumption fails for thecases where the media contained by the coarse cells is highly heterogeneous,as it is well known that heterogeneous conductive media induce a non-linear andnon-smooth behavior of the electric field across media interfaces ([169]). In par-ticular, when the heterogeneity is located close to the boundary of the coarse cell,the non-linear behavior of the electric field significantly violates the assumptionof linear fields at the boundaries. Hence, by imposing linear boundary conditionsin such cases for the construction of the multiscale basis functions, the MSFVmethod creates a mismatch between the true and the multiscale solution acrossthe coarse cell boundary. This mismatch yields to produce inaccurate solutions.In this section, we propose an oversampling technique to overcome this difficulty.Oversampling methods are used to reduce boundary effects in the construc-tion of the multiscale basis functions per single coarse-mesh cell [51, 89]. Themain idea is to compute the multiscale basis functions using a local extendeddomain, and to use only the fine-mesh information at the interior of the cell toconstruct the multiscale basis functions.We now proceed to develop our oversampling technique. To do so, we adaptthe oversampling technique originally proposed by [89] for linear elliptic problemswith strongly discontinuous coefficients and for a nodal FE discretization, to ap-ply to the MSFV method for EM modeling with edge variables discussed in theprevious section, which uses a staggered FV discretization (i.e., the MFV methoddiscussed in Section 2.4).For a given coarse cell ΩHk , the core idea behind our oversampling methodconsists of the following two steps, which are illustrated in Figure 4.2.First, we compute multiscale basis functions using a local extended domain,denoted as ΩH,extk , which includes the coarse cell ΩHk and a neighborhood of fine88Figure 4.2: Schematic representation of the two main steps to implementthe oversampling method.cells around it. If the coarse cell ΩHk is at the boundary of the computationaldomain, then we only extend the local domain to where it is possible. To computethe multiscale basis functions in ΩH,extk , we formulate the twelve local Maxwellsystems as in (4.1)-(4.3), but rather than using ΩHk as the local domain we useΩH,extk , then we apply the MFV method as discussed in Section 2.4. We referto the set of discrete solutions {ek,ext1 , . . . ,ek,ext12 } as extended multiscale basisfunctions for the cell ΩHk .Second, we use the set of extended multiscale basis functions obtained in theprevious step to compute the actual set of multiscale basis functions {ek1, . . . ,ek12}in ΩHk . Since the construction of the multiscale basis is done cell by cell, thereis no guaranty that the tangential components of the multiscale basis functionsare continuous at the boundary of the coarse cell ΩHk . In order to mitigate thisissue, we impose the following weak-continuity condition in the construction of the89multiscale basis functions to guarantee they will be weakly continuous along eachshared boundary among immediate neighboring coarse cells,Aedgem(~Ekl):=1Ledgem∫edgem~Ekl ·~τedgem ds = δml; m, l = 1, . . . ,12, (4.8)where ~Ekl denotes the continuous form of the lth multiscale basis function ekl ,Ledgem denotes the length of the mth edge of ΩHk , ~τedgem denotes the unit tan-gent vector to the mth edge of ΩHk , and δml is the Kronecker delta. That is, wetake a ‘normalized average’ of the multiscale basis functions at the boundary ofthe coarse cell. This condition is equivalent to the definition of edge degrees offreedom of a staggered cell in the context of Ne´de´lec finite elements [130, 131].Note the difference with the continuity condition (4.4) imposed in the constructionof the multiscale functions of the MSFV method without oversampling. Integratingnumerically the continuity condition (4.8), we can express it asˆAedgem(ekl)≈ v>edgemekl = δml; m, l = 1, . . . ,12, (4.9)where vedgem is the vector that computes the normalized line integral along themth edge of ΩHk .Using (4.9) and following the main lines given in the oversampling techniqueproposed by [89], we continue the development of our oversampling technique byshowing how to compute {ek1, . . . ,ek12} from {ek,ext1 , . . . ,ek,ext12 } in detail.We begin by expressing the jth multiscale basis function, ekj, as a linear com-bination of the set of extended basis functions as followsekj =12∑l=1cl, j ek,extl = [ek,ext1 , . . . ,ek,ext12 ]c j; j = 1, . . . ,12, (4.10)where c j = [c1, j, . . . ,c12, j]> are coefficients to be determined. Now, to determineuniquely such coefficients, we apply condition (4.9) to (4.10), which results in the90following system of equationsˆAedge1(ek,ext1 ) . . .ˆAedge1(ek,ext12 )ˆAedge2(ek,ext1 ) . . .ˆAedge2(ek,ext12 ).... . ....ˆAedge12(ek,ext1 ) . . .ˆAedge12(ek,ext12 )c1,1 . . . c1,12c2,1 . . . c2,12.... . ....c12,1 . . . c12,12= I12×12,(4.11)where I12×12 denotes the 12×12 identity matrix. Combining equations (4.9),(4.10) and (4.11), we obtain the expression for the desired coefficients, that isC=v>edge1ek,ext1 . . . v>edge1ek,ext12v>edge2ek,ext1 . . . v>edge2ek,ext12.... . ....v>edge12ek,ext1 . . . v>edge12ek,ext12−1. (4.12)Now, C is invertible because its columns are linearly independent. To see this,we consider the matrix (4.12) before it is inverted. To show that the columns ofthis matrix form a linearly independent set, we need to show that for any linearcombination of such vectors:α1v>edge1ek,ext1v>edge2ek,ext1...v>edge12ek,ext1+ · · ·+α12v>edge1ek,ext12v>edge2ek,ext12...v>edge12ek,ext12=00...0 (4.13)the scalars α1, . . . ,α12 are zero. As described before (4.9), v>edgem ; m = 1, . . . ,12denotes the vector that computes the normalized line integral along the mth edgeof the coarse cell ΩHk . Now, we can write expression (4.13) asα1v>edge1v>edge2...v>edge12ek,ext1 + · · ·+α12v>edge1v>edge2...v>edge12ek,ext12 =v>edge1v>edge2...v>edge1212∑l=1αlek,extl =~0.(4.14)91Since {ek,extl }12l=1 form a basis, then αl = 0, ∀l = 1, . . . ,12.After we construct the multiscale basis functions {ek1, . . . ,ek12} using our over-sampling technique, we continue to follow the procedure for the MSFV method(Figure 4.1) to compute the solution. That is, this set of twelve multiscale basisfunctions enable the use of the local interpolation matrix Pk, given by (4.5), withinthe assembly of the global coarse-to-fine interpolation matrix P. The interpola-tion matrix P is then used within a Galerkin formulation to obtain the coarse-meshsystem (4.6), which we ultimately solve.4.4 Numerical Results in 3DIn this section, we demonstrate the accuracy and computational performance ofour multiscale finite volume method with oversampling (MSFV+O) by simulatingEM responses for two 3D synthetic electrical conductivity models: one with amineral deposit in a geologically-rich medium and one with a random isotropicheterogeneous medium. As analytical solutions are not available for these twoexamples, the results of these simulations are compared to the simulation resultsfrom the fine-mesh reference models, respectively. Since we need to use a finemesh to compute the multiscale basis functions numerically, it makes sense tocompare our MSFV+O method with a traditional FV method at the fine mesh. Thisexample also demonstrates the power of combining adaptive mesh refinementtechniques with a multiscale approach to produce faster and accurate simulations.4.4.1 Simulations for the Synthetic Lalor ModelFor the first example, we use the synthetic Lalor electrical conductivity modelintroduced in Section 3.5.2, which is based on the inversion results of field mea-surements over the Canadian Lalor mine obtained by [177].As discussed in Section 3.5.2, the Lalor model, shown in Figure 4.3, hasan area with non-flat topography and extends from 0 to 6.5 km along the x, yand z directions, respectively. The model comprises air and the subsurface thatis composed of 35 geologic units. The unit with the largest conductivity valuerepresents the mineral deposit, which is composed of three bodies. We assumea conductivity of 10−8 S/m in the air. The subsurface conductivity values range92Figure 4.3: Subsurface part of the synthetic electrical conductivity Lalormodel and large-loop EM survey setup. The model is discretized ona fine OcTree mesh with 546,295 cells. The conductivity varies overfive orders of magnitude throughout the whole model.from 1.96× 10−5 to 0.28 S/m. The magnetic permeability takes its free spacevalue (i.e., µ = 4pi×10−7 Vs/Am).We consider a large-loop EM survey for this example, where we use a rectan-gular transmitter loop with dimensions 2 km × 3 km, operating at the frequenciesof 1, 10, 20, 40, 100, 200 and 400 Hz. The transmitter is placed on the earth’ssurface and it is centered above the largest body of the mineral deposit, as shownin Figure 4.3. Inside the loop, we place a uniform grid of receivers that measurethe three components of the magnetic flux (~B = [Bx,By,Bz]>). The receivers areseparated by 50 m along the x and y directions, respectively.Our aim is to estimate the secondary magnetic flux induced by the mineraldeposit in the survey area. For this purpose, we simulate two sets of the magnetic93flux data for each frequency. The first data set considers the conductivity modelincluding all geologic units, and the second data set excludes the mineral depositfrom the original conductivity model. Each of these two data sets consists of themeasurements of ~B taken at the receiver locations. The secondary magnetic fluxinduced by the mineral deposit at the survey area, denoted as ∆~Bdeposit, is thencomputed by subtracting the two data sets. The estimation of ∆~Bdeposit requirescarrying out two forward model simulations per frequency.To compute a reference solution, we discretize the conductivity model usingthe MFV method as discussed in [72, 86] at the same fine OcTree mesh describedin Section 3.5.2, which is shown in Figure 4.3. The fine mesh has 546,295 cells.Considering the averaged background conductivity value (4.5 ×10−3 S/m), theskin depths for the frequencies of 1, 10, 20, 40, 100, 200 and 400 Hz are roughly7,498, 2,371, 1,677, 1,186, 750, 530 and 375 m, respectively. Thus, using cellssizes of (50 m)3 within the survey area continues to be sufficient to capture thebehavior of the EM fields in this setup. Using the MFV method on the fine OcTreemesh yields systems with roughly 1.5 millions edge DOF which we solve usingthe parallel sparse direct solver MUMPS [3]. The average computation time persingle simulation is roughly 721 s on a two hexa-core Intel Xeon X5660 CPUs at2.8 Hz with 64 GB shared RAM using the parallel computing toolbox of MATLAB([162]). Figure 4.4 shows the Euclidean norm of the total, real and imaginaryparts of ∆~Bdeposit for each frequency considered. The real and imaginary parts ofthe results obtained for the z-component of ∆~Bdeposit, denoted as ∆Bzdeposit, at 100 Hzare shown in Figure 4.6(a) and Figure 4.7(a), respectively.In order to use the MSFV+O method introduced in Section 4.3, we need tochoose a suitable coarse mesh to discretize the conductivity model and the sizeof the local extended domain to compute its corresponding projection matrix.As a coarse mesh, we consider the same coarse OcTree mesh nested in thefine OcTree mesh previously described in Section 3.5.2, which is shown in Figure4.5. The coarse OcTree mesh is designed to maintain the fine-mesh resolutionof (50 m)3 inside the survey area, whereas the rest of the domain is filled withgradually increasing coarser cells. In total, this coarse mesh has 60,656 cells. Toanalyze the performance of our MSFV+O method for coarse OcTree meshes, wedo not refine the mesh outside the survey area where a large conductivity contrast94Figure 4.4: Euclidean norm of total, real and imaginary parts of the ref-erence (fine-mesh) solution for the Lalor conductivity model per fre-quency.is present in the model. For example, this mesh discretizes the mineral depositwith cells of size (200 m)3 and (400 m)3, and the non-flat topography with cellsof size (400 m)3 and (800 m)3. We expect the simulations for the frequenciesof 200 and 400 Hz to be particularly challenging for our MSFV+O method as thecoarsening in those areas can be considered extreme due to the cell size is in theorder of the skin depth.Next, to investigate the effect of the size of the local extended domain, i.e., thenumber of fine-mesh padding cells by which we extend every coarse cell, on theresulting magnetic flux data, we pad the coarse cell using 2, 4 and 8 fine cells.Each fine padding cell is of size (50 m)3. The chosen local extended domain sizescorrespond to extending each (200 m)3-coarse cell by half, one and two coarsecells, respectively. The (200 m)3-coarse cells are the majority of the coarse cellswhere the largest conductivity contrast happen in this setting (Figure 4.5).Applying MSFV and MSFV+O on the coarse mesh shown in Figure 4.5, we95Figure 4.5: Subsurface part of the synthetic electrical conductivity Lalormodel and large-loop EM survey setup. The model is discretized ona coarse OcTree mesh with 60,656 cells. The conductivity varies overeight orders of magnitude throughout the whole model.obtain reduced linear systems with 169,892 DOF, which are also solved usingMUMPS. When using MSFV+O the total average run times per single simulationfor extended domain sizes of 2, 4 and 8 padding cells are roughly 185 s, 442s and 1 h, respectively, on the same machine. The real and imaginary partsof ∆Bzdeposit at 100 Hz for extended domain sizes of 2, 4 and 8 padding cells areshown in Figures 4.6(b), 4.6(c) and 4.6(d), respectively; and Figures 4.7(b), 4.7(c)and 4.7(d), respectively. In order to use MSFV (without oversampling), we firstadapt this method for OcTree meshes, as the original version is derived for tensormeshes only (cf. [73]). In this case, the total average run time per single simulationis roughly 73 s on the same machine. The real and imaginary parts of ∆Bzdeposit at100 Hz are shown in Figure 4.6(e) and Figure 4.7(e), respectively.96We also carry out MFV simulations using homogenized electrical conductivitymodels that we construct using volume-arithmetic, -geometric and -harmonic av-eraging of the fine-mesh conductivity inside each coarse cell of the OcTree meshshown in Figure 4.5. Doing so allows us to compare the accuracy that MFV so-lutions achieve on the coarse mesh with the one achieved by MSFV with andwithout oversampling. The total average run time per single simulation is roughly128 s on the same machine. The real and imaginary parts of ∆Bzdeposit at 100 Hz foreach of the three homogenized solutions are shown in Figures 4.6(f), 4.6(g) and4.6(h), respectively; and Figures 4.7(f), 4.7(g) and 4.7(h), respectively.Figures 4.6 and 4.7 show the results obtained for the real and imaginary partsof ∆Bzdeposit with MFV on the fine-mesh, MSFV, MSFV+O with three different over-sampling sizes, and MFV with three different homogenized conductivity models forthe frequency of 100 Hz, respectively. All results are plotted using the same colorscale and range. From these figures we see that our oversampling technique pro-duces the most accurate results in comparison with the rest of the methods used.Observe that using only an oversampling size of 2 padding cells significantly im-proves the quality of the approximations, whereas MSFV tends to overestimatethe approximation for both real and imaginary parts of ∆Bzdeposit.97(a) (b) (c) (d)(e) (f) (g) (h)Figure 4.6: Real part of the z component of the secondary magnetic flux induced by the mineral deposit in the surveyarea, ∆Bzdeposit, for our large-loop EM survey at 100 Hz. (a) reference solution computed using the MFV methodon the fine OcTree mesh with 546,295 cells. (b), (c) and (d): results using the MSFV+O method with 2, 4 and8 padding cells on the coarse OcTree mesh, respectively. (e): results using the MSFV (without oversampling)method on the coarse OcTree mesh with 60,656 cells. (f), (g) and (h): results using MFV with the conduc-tivity model homogenized using arithmetic, geometric and harmonic averaging on the coarse OcTree mesh,respectively. All results are shown in picoteslas (pT) and plotted using the same color scale.98(a) (b) (c) (d)(e) (f) (g) (h)Figure 4.7: Imaginary part of the z component of the secondary magnetic flux induced by the mineral deposit inthe survey area, ∆Bzdeposit, for our large-loop EM survey at 100 Hz. (a) reference solution computed using theMFV method on the fine OcTree mesh with 546,295 cells. (b), (c) and (d): results using the MSFV+O methodwith 2, 4 and 8 padding cells on the coarse OcTree mesh, respectively. (e): results using the MSFV (withoutoversampling) method on the coarse OcTree mesh with 60,656 cells. (f), (g) and (h): results using MFV withthe conductivity model homogenized using arithmetic, geometric and harmonic averaging on the coarse OcTreemesh, respectively. All results are shown in picoteslas (pT) and plotted using the same color scale.99Table 4.1 shows the relative errors in Euclidean norm for the total, real andimaginary parts of ∆~Bdeposit obtained from comparing the reference (fine-mesh)solution with the MSFV, MSFV+O and MFV with three different homogenized so-lutions for each frequency and local extended domain size. Table 4.2 shows therelative errors in Euclidean norm for the magnitude of the three components of∆~Bdeposit (∆~Bdeposit = [∆Bxdeposit,∆Bydeposit,∆Bzdeposit]>) obtained from comparing the refer-ence (fine-mesh) solution with the MSFV, MSFV+O and MFV with three differenthomogenized solutions for each frequency and local extended domain size. Therelative error is computed as the ratio of the Euclidean norm of the difference ofthe fine and coarse-mesh data to the Euclidean norm of the fine-mesh data. Thecoarse-mesh solutions are not interpolated to the fine mesh to compute the error.Table 4.3 summarizes the average run time per single simulation required to com-pute ∆~Bdeposit for each of the methods discussed. From these tables we observethe following:First, we see that our oversampling technique significantly improves the ac-curacy for the total, real and imaginary response as well as for each of the threecomponents of the solution in comparison to the MSFV and MFV with three differ-ent homogenized conductivity models as the errors decrease with oversampling.In particular, it is surprising to see how well MFV with simple geometric averagingdid when compared with MSFV.Second, as the size of the local extended domain increases the error de-creases at the expense of more computational run time, which, however, is stillconsiderably lower compared to the time of the reference solution for the cases of2 and 4 padding cells (see Table 4.3). These results suggest that by using a localextended domain size of at least a half the number of fine cells contained in thecoarse cell(s) where the major contrast of conductivity happens, we may increasesignificantly the accuracy obtained with MSFV+O.Third, the magnitude of the errors come from comparing secondary magneticfield data obtained on the fine and coarse mesh, respectively. The data is theresponse of the deep and buried conductive mineral deposit, which constitutes arather weak secondary field compared to the primary field induced by the air andbackground conductivity. The level of accuracy of this simulation is relative to thetotal signal, which values are quite small (see Figure 4.4). We can improve the100accuracy of MSFV and MSFV+O by refining the coarse mesh and interpolating thecoarse-mesh solution to the fine mesh as shown in [73]. Yet, MSFV+O providesa reasonable approximation to the fine-mesh solution using a coarse mesh that isonly 10% the size of the fine mesh.Fourth, the errors for the imaginary part of ∆~Bdeposit for the frequencies of 1, 10,20 and 40 Hz resemble the errors for total ∆~Bdeposit due to the real part of ∆~Bdepositis up to two orders of magnitude smaller than their corresponding imaginary parts(Figure 4.4). For the frequencies of 100 and 400 Hz, where the real and imaginaryparts of ∆~Bdeposit are roughly in the same order of magnitude, we see the relativeerror for total ∆~Bdeposit represents the error contributions of these two componentsof the data.Fifth, for the simulation at 200 Hz the errors for the real part of ∆~Bdeposit resem-ble the errors for total ∆~Bdeposit due to the imaginary part of ∆~Bdeposit is one order ofmagnitude smaller (Figure 4.4). The large errors in the imaginary part of the dataas well as the increase and decrease in the error going from padding cell 2 to 4and 4 to 8, respectively, can be attributed to the combined discretization error insimulating secondary field data and the extreme coarsening in the mesh shown inFigure 4.5. For this frequency, the cell sizes used to discretize the mineral depositare in the order of the skin depth.Sixth, for a simulation at 400 Hz with an extended domain size of 8 paddingcells the slight increment in the error of the real part may be attributed to the ex-cessive coarsening in the mesh for such high frequency. In this case, the coarsecells used to discretize the mineral deposit are larger than the skin depth. Despitethe excessive coarsening, our oversampled approximation yields comparable re-sults to the reference (fine-mesh) solution for the largest two frequencies.Seventh, the increment of the MSFV+O setup time as the size of the localextended domain increases (Table 4.3) comes from constructing the local inter-polation matrices. Since the computation of these matrices is done locally insideeach coarse cell independently of each other, this process can be done in par-allel (see discussion in Section 4.2). A robust implementation using a parallelcommunication protocol (e.g. Message Passage Interface (MPI)) for program-ming parallel computers should greatly reduce the overhead time in constructingthese interpolation matrices. The research work in [46], which focuses in prob-101lems in porous media, shows progress along this lines. However, we note thatusing oversampling with 2 padding cells is already increasing the accuracy of thesolution and the computational time is significantly lower that that of computingthe reference solution.Table 4.1: Relative errors in Euclidean norm for the total, real and imaginaryparts of the secondary magnetic fluxes induced by the mineral depositin the survey area, ∆~Bdeposit. Note that pc stands for padding cells.Table of relative errors in Euclidean norm for ∆~BdepositFrequencyMethod 1Hz 10Hz 20Hz 40Hz 100Hz 200Hz 400HzRelative errors for total ∆~Bdeposit (per cent)MFV+Arithmetic 192.05 191.60 190.43 187.70 183.25 171.94 149.34MFV+Geometric 68.57 68.51 68.34 67.72 64.96 60.61 55.62MFV+Harmonic 97.20 97.20 97.18 97.11 96.81 96.30 95.57MSFV 69.70 69.72 69.79 70.32 73.09 72.65 63.93MSFV+O (2 pc) 15.84 15.83 15.79 15.75 16.17 15.98 13.46MSFV+O (4 pc) 13.31 13.33 13.38 13.60 14.48 14.63 11.36MSFV+O (8 pc) 10.63 10.65 10.72 10.98 12.20 12.67 10.08Relative errors for real part of ∆~Bdeposit (per cent)MFV+Arithmetic 216.47 214.97 211.12 202.06 188.63 170.99 159.57MFV+Geometric 81.01 80.95 80.75 79.72 73.05 60.33 34.18MFV+Harmonic 98.46 98.45 98.43 98.33 97.74 96.52 93.45MSFV 73.09 72.92 72.56 72.29 74.22 70.82 59.47MSFV+O (2 pc) 21.41 21.34 21.16 20.57 18.11 14.21 8.25MSFV+O (4 pc) 18.36 18.38 18.39 18.19 15.92 12.49 8.14MSFV+O (8 pc) 16.12 16.10 16.02 15.49 12.74 10.40 8.51Relative errors for imaginary part of ∆~Bdeposit (per cent)MFV+Arithmetic 192.05 191.26 189.26 184.63 174.80 197.75 133.27MFV+Geometric 68.57 68.32 67.61 64.99 50.30 68.40 76.55MFV+Harmonic 97.20 97.18 97.11 96.86 95.38 89.47 98.52MSFV 69.70 69.67 69.64 69.91 71.34 114.01 69.83MSFV+O (2 pc) 15.84 15.74 15.45 14.57 12.68 42.77 18.53MSFV+O (4 pc) 13.31 13.24 13.06 12.46 11.99 43.84 14.81MSFV+O (8 pc) 10.63 10.55 10.36 9.81 11.35 41.16 11.98102Table 4.2: Relative errors (per cent) in Euclidean norm for the magnitude ofthe x, y and z components of the secondary magnetic fluxes induced bythe mineral deposit in the survey area, ∆~Bdeposit. Note that pc stands forpadding cells.Method ∆Bxdeposit ∆Bydeposit ∆Bzdeposit ∆Bxdeposit ∆Bydeposit ∆BzdepositFrequency of 1 HzMFV+Arithmetic 176.44 202.20 195.12MFV+Geometric 66.26 71.52 68.65MFV+Harmonic 97.03 97.38 97.22MSFV 65.91 70.44 70.95MSFV+O (2 pc) 13.37 19.83 15.53MSFV+O (4 pc) 12.04 16.24 12.90MSFV+O (8 pc) 10.45 12.13 10.25Frequency of 10 Hz Frequency of 100 HzMFV+Arithmetic 176.20 201.86 194.59 170.48 198.30 184.62MFV+Geometric 66.17 71.50 68.60 62.34 68.63 65.11MFV+Harmonic 97.03 97.37 97.21 96.63 97.05 96.82MSFV 66.06 70.39 70.94 70.39 73.56 74.16MSFV+O (2 pc) 13.39 19.78 15.51 14.00 19.41 16.13MSFV+O (4 pc) 12.05 16.25 12.93 12.56 17.12 14.52MSFV+O (8 pc) 10.45 12.14 10.29 11.26 13.31 12.30Frequency of 20 Hz Frequency of 200 HzMFV+Arithmetic 175.59 201.03 193.21 159.79 189.57 172.55MFV+Geometric 65.91 71.44 68.44 58.65 63.69 60.68MFV+Harmonic 97.00 97.37 97.20 96.20 96.47 96.30MSFV 66.48 70.30 70.95 68.82 75.18 73.73MSFV+O (2 pc) 13.45 19.65 15.48 13.41 19.21 16.15MSFV+O (4 pc) 12.07 16.28 13.01 12.76 16.76 14.85MSFV+O (8 pc) 10.47 12.18 10.39 10.48 15.03 12.92Frequency of 40 Hz Frequency of 400 HzMFV+Arithmetic 174.12 199.32 189.94 141.64 163.62 148.96MFV+Geometric 65.09 71.10 67.86 54.59 57.49 55.60MFV+Harmonic 96.92 97.33 97.13 95.57 95.60 95.55MSFV 67.73 70.43 71.36 60.31 67.47 64.68MSFV+O (2 pc) 13.63 19.38 15.48 11.04 16.43 13.64MSFV+O (4 pc) 12.15 16.42 13.33 9.52 13.52 11.56MSFV+O (8 pc) 10.56 12.36 10.76 7.83 12.82 10.23103Table 4.3: Average run time in seconds per simulation required to compute∆~Bdeposit on a two hexa-core Intel Xeon X5660 CPUs at 2.8 Hz with 64GB shared RAM using MATLAB. Setup time: time required to computethe local interpolation matrices and to assemble the reduced system ofequations to be solved. Solve time: time to solve the reduced system ofequations. Total time: the sum of the setup and solve times.Method Setup time (s) Solve time (s) Total time (s)MFV on fine OcTree mesh - 721 s 721 sMFV + Arithmetic - 127 s 127 sMFV + Geometric - 128 s 128 sMFV + Harmonic - 128 s 128 sMSFV 42 s 31 s 73 sMSFV+O (2 padding cells) 147 s 38 s 185 sMSFV+O (4 padding cells) 410 s 32 s 442 sMSFV+O (8 padding cells) 3,566 s 32 s 3,598 sThe next example demonstrates the effect of considering a different heteroge-neous conductivity model for the same survey and meshes setup on the resultsproduced by our proposed MSFV+O method.4.4.2 Simulations for Random Heterogeneous Isotropic MediaIn this section, we construct an electrical conductivity model of a random isotropicheterogeneous medium to provide a second magnetic flux data set to validate theproposed oversampling technique. For this example, we use the same large-loopEM survey configuration, computational domain and fine and coarse-mesh setupas described in Section 4.4.1.The new synthetic conductivity model, shown in Figure 4.8, also comprises airand the subsurface that is composed of random heterogeneous isotropic media.We assume that the subsurface conductivity is log(10)-normally distributed withmean of 2.5× 10−3 and standard variation of 0.4. To reuse the fine-mesh setupdescribed in Section 4.4.1, we mapped the values of the subsurface conductivityto the interval [10−5,10−1] S/m. Considering the averaged background conduc-tivity value (1.6 ×10−3 S/m), the skin depths for the frequencies of 1, 10, 20, 40,100, 200 and 400 Hz are roughly 12,575, 3,977, 2,812, 1,988, 1,258, 889 and104629 m, respectively. Thus, using cells sizes of (50 m)3 within the survey areacontinues to be sufficient to capture the behavior of the EM fields in this setup.As in the previous example, we assume that the conductivity of the air is of 10−8S/m and that the magnetic permeability takes its free space value. Note that thisconductivity model does not contain clearly defined conductive features (as in theLalor example); instead the contrast of conductivity is distributed throughout theentire subsurface part of the model. This allows us to challenge the oversamplingtechnique with large conductivity contrast across the material interfaces.The goal of the simulation in this case is to compute the secondary magneticflux induced by the random heterogeneous conductive medium in the survey area.To do so, we simulate two sets of the magnetic flux data for each frequency. Thefirst data set considers the conductivity model including air and the subsurface;the second data set considers a conductivity model where the subsurface is re-placed by air. Each of these two data sets consists of the measurements of thethree components of ~B taken at the receiver locations ([~B = Bx,By,Bz]>). Thesecondary magnetic flux induced by the random conductive medium at the surveyarea, denoted as ∆~Brandom, is then computed by subtracting the two data sets. Theestimation of ∆~Brandom requires two forward model simulations per frequency.Now, to complete the validation for our oversampling technique we need tocompute a reference solution, a multiscaled solution with and without oversam-pling and averaged-based homogenized solutions. We use the same machinedescribed in Section 4.4.1, MUMPS and MATLAB to run all the simulations pre-sented in this section.To compute a reference solution, we apply the MFV method on the staggeredfine OcTree mesh shown in Figure 4.8 (see [72, 86] for details). The averagecomputation time per single simulation is roughly 712 s. Figure 4.9 shows theEuclidean norm of the total, real and imaginary parts of ∆~Brandom for each fre-quency considered. The real and imaginary parts of the results obtained for the z-component of ∆~Brandom, denoted as ∆Bzrandom, at 100 Hz are shown in Figures 4.11(a)and 4.12(a), respectively.Next, we apply the MSFV method with and without oversampling on the coarseOcTree mesh shown in Figure 4.10 to simulate ∆~Brandom. For this example, we useextended domain sizes of 1, 2 and 4 fine-mesh padding cells. Each fine padding105Figure 4.8: Subsurface part of the random heterogeneous isotropic electri-cal conductivity model and large-loop EM survey setup. The model isdiscretized on a fine OcTree mesh with 546,295 cells. The conductivityvaries over eight orders of magnitude throughout the whole model.cell is of size (50 m)3. When using MSFV+O the total average run times per singlesimulation for extended domain sizes of 1, 2 and 4 padding cells are roughly 106,156 and 439 s, respectively. The real and imaginary parts of ∆Bzrandom at 100 Hz forextended domain sizes of 1, 2 and 4 padding cells are shown in Figures 4.11(b),4.11(c) and 4.11(d), respectively; and Figures 4.12(b), 4.12(c) and 4.12(d), re-spectively. On the other hand, when using MSFV (without oversampling), the totalaverage run time per single simulation is roughly 74 s. The real and imaginary106Figure 4.9: Euclidean norm of total, real and imaginary parts of the refer-ence (fine-mesh) solution for the random conductivity model per fre-quency.parts of ∆Bzrandom at 100 Hz are shown in Figures 4.11(e) and 4.12(e), respectively.Finally, we carry out MFV simulations using homogenized electrical conductiv-ity models that we construct using volume-arithmetic, -geometric and -harmonicaveraging of the fine-mesh conductivity inside each coarse cell of the OcTreemesh shown in Figure 4.10. The average run time per single simulation is roughly115 s. The real and imaginary parts of ∆Bzrandom at 100 Hz for each of the threehomogenized solutions are shown in Figures 4.11(f), 4.11(g) and 4.11(h), respec-tively; and Figures 4.12(f), 4.12(g) and 4.12(h), respectively.Table 4.4 shows the relative errors in Euclidean norm for the total, real andimaginary parts of ∆~Brandom obtained from comparing the reference (fine-mesh)solution with the MSFV, MSFV+O and MFV with three different homogenized so-lutions for each frequency and local extended domain size. Table 4.5 shows therelative errors in Euclidean norm for the magnitude of the three components of∆~Brandom (∆~Brandom = [∆Bxrandom,∆Byrandom,∆Bzrandom]>) obtained from comparing the ref-107Figure 4.10: Subsurface part of the random heterogeneous isotropic electri-cal conductivity model and large-loop EM survey setup. The modelis discretized on a coarse OcTree mesh with 60,656 cells. The con-ductivity varies over eight orders of magnitude throughout the wholemodel.erence (fine-mesh) solution with the MSFV, MSFV+O and MFV with three differenthomogenized solutions for each frequency and local extended domain size. Onceagain, the relative error is computed as the ratio of the Euclidean norm of the dif-ference of the fine and coarse-mesh data to the Euclidean norm of the fine-meshdata. The coarse-mesh solutions are not interpolated to the fine mesh to com-pute the error. Table 4.6 summarizes the average run time per single simulationrequired to compute ∆~Brandom for each of the methods discussed. Figures 4.11 and1084.12 show the results obtained for the real and imaginary parts of ∆Bzrandom withMFV on the fine-mesh, MSFV, MSFV+O with three different oversampling sizes,and MFV with three different homogenized conductivity models for the frequencyof 100 Hz, respectively. All results are plotted using the same color scale andrange. From these tables and these figures we observe the following:First, once again, we see that our oversampling technique improves the ac-curacy for the total, real and imaginary response as well as for each of the threecomponents of the solution in comparison to the MSFV and MFV with three dif-ferent homogenized conductivity models as the errors decrease with oversam-pling. As expected, the homogenized conductivity model obtained from usingthe volume-geometric mean produces the most accurate results as compared tothose obtained using volume-arithmetic and volume-harmonic means. The geo-metric mean provides a reliable approximation of the effective conductivity if thefine-scale variation satisfies certain conditions [163]. We also see that the multi-scaled solutions are more accurate than the homogenized solutions.Second, using a local extended domain size of 1 fine padding cell suffices toimprove the accuracy of the solution when compared with the rest of the methodsdiscussed and it is slightly more expensive than the solution obtained with MSFV.We see that as the size of the local extended domain increases the error alsowas slightly increased. This small increment in the error is still lower comparedto the error of the MSFV solution for most of the cases, except for the imaginarypart of ∆~Brandom at 100 Hz and the real part of ∆~Brandom at 400 Hz. We attribute theincrement in the error to the random and uncorrelated structure of the conductivitymodel. On the one hand, increasing the oversampling size reduces the effect ofthe imposed linear boundary conditions. This explains the improved accuracy ofthe MSFV+O over MSFV. On the other hand, increasing the oversampling sizetakes into account structures outside the core coarse cell that are not necessarilycorrelated to the structures inside the core coarse cell. This biases the multiscalebasis towards structures outside the core cell and decreases the local approxima-tion quality.Third, comparing the relative accuracy for the two examples presented wenote the two data types differ significantly. In the Lalor example (Section 4.4.1),the data is the response of the deep and buried conductive mineral deposit, which109constitutes a rather weak secondary field compared to the primary field induced bythe air and background conductivity (Figure 4.4). In the random medium example,the data is the response of the entire subsurface; the primary field here is onlyfor the air (Figure 4.9). Therefore, the ratio of the secondary to primary field ismuch higher in the random example than for the Lalor example. Since the levelof accuracy of a simulation is relative to the total signal, the accuracy of a weakersecondary signal is lower than for a stronger secondary signal.Fourth, as in the previous example, the errors for the imaginary part of ∆~Brandomfor the frequencies of 1, 10, 20 and 40 Hz resemble the errors for total ∆~Brandom dueto the real part of ∆~Brandom is considerably smaller than their corresponding imag-inary parts (Figure 4.9). For the frequencies of 100, 200 and 400 Hz, where thereal and imaginary parts of ∆~Brandom are roughly in the same order of magnitude(Figure 4.9), we see the relative error for total ∆~Brandom represents the error contri-butions of these two parts of the data.110(a) (b) (c) (d)(e) (f) (g) (h)Figure 4.11: Real part of the z component of the secondary magnetic flux induced by the random heterogeneousconductive medium in the survey area, ∆Bzrandom, for our large-loop EM survey at 100 Hz. (a) reference solutioncomputed using the MFV method on the fine OcTree mesh with 546,295 cells. (b), (c) and (d): results usingthe MSFV+O method with 1, 2 and 4 padding cells on the coarse OcTree mesh, respectively. (e): results usingthe MSFV (without oversampling) method on the coarse OcTree mesh with 60,656 cells. (f), (g) and (h): resultsusing MFV with the conductivity model homogenized using arithmetic, geometric and harmonic averaging onthe coarse OcTree mesh, respectively. All results are shown in picoteslas (pT) and plotted using the same colorscale.111(a) (b) (c) (d)(e) (f) (g) (h)Figure 4.12: Imaginary part of the z component of the secondary magnetic flux induced by the random heterogeneousconductive medium in the survey area, ∆Bzrandom, for our large-loop EM survey at 100 Hz. (a) reference solutioncomputed using the MFV method on the fine OcTree mesh with 546,295 cells. (b), (c) and (d): results usingthe MSFV+O method with 1, 2 and 4 padding cells on the coarse OcTree mesh, respectively. (e): results usingthe MSFV (without oversampling) method on the coarse OcTree mesh with 60,656 cells. (f), (g) and (h): resultsusing MFV with the conductivity model homogenized using arithmetic, geometric and harmonic averaging onthe coarse OcTree mesh, respectively. All results are shown in picoteslas (pT) and plotted using the same colorscale.112Table 4.4: Relative errors in Euclidean norm for the total, real and imaginaryparts of the secondary magnetic fluxes induced by the random hetero-geneous conductive medium in the survey area, ∆~Brandom. Note that pcstands for padding cells.Table of relative errors in Euclidean norm for ∆~BrandomFrequencyMethod 1 Hz 10 Hz 20 Hz 40 Hz 100 Hz 200 Hz 400 HzRelative errors for total ∆~Brandom (per cent)MFV + Arithmetic 10.50 9.95 8.67 6.27 3.90 2.69 1.81MFV + Geometric 9.77 9.51 8.83 7.24 4.96 3.35 1.96MFV + Harmonic 21.38 20.94 19.76 16.53 10.64 7.09 3.84MSFV 4.41 4.22 3.74 2.77 1.75 1.20 0.76MSFV+O (1 pc) 0.43 0.53 0.66 0.80 0.66 0.54 0.44MSFV+O (2 pc) 0.72 0.80 0.94 1.05 0.83 0.64 0.49MSFV+O (4 pc) 0.93 0.99 1.10 1.17 0.87 0.66 0.49Relative errors for real part of ∆~Brandom (per cent)MFV + Arithmetic 38.36 35.03 27.89 16.56 7.68 3.52 0.75MFV + Geometric 25.39 24.38 21.88 16.32 9.33 5.32 2.53MFV + Harmonic 53.99 52.57 48.75 38.45 20.64 11.03 4.02MSFV 14.96 13.89 11.46 7.14 3.44 1.67 0.44MSFV+O (1 pc) 1.60 1.32 0.71 0.72 1.00 0.78 0.59MSFV+O (2 pc) 0.71 0.50 0.50 1.37 1.38 0.97 0.65MSFV+O (4 pc) 0.44 0.57 1.02 1.79 1.51 1.02 0.67Relative errors for imaginary part of ∆~Brandom (per cent)MFV + Arithmetic 10.49 9.11 6.32 2.59 0.44 1.99 2.51MFV + Geometric 9.77 9.12 7.62 4.79 1.78 0.78 1.00MFV + Harmonic 21.37 20.14 17.06 10.31 2.35 2.39 3.63MSFV 4.41 3.91 2.85 1.29 0.24 0.79 0.99MSFV+O (1 pc) 0.43 0.51 0.67 0.81 0.49 0.31 0.19MSFV+O (2 pc) 0.73 0.81 0.95 1.00 0.52 0.28 0.19MSFV+O (4 pc) 0.93 1.00 1.10 1.04 0.48 0.25 0.20113Table 4.5: Relative errors (per cent) in Euclidean norm for the magnitude ofthe x, y and z components of the secondary magnetic fluxes induced bythe mineral deposit in the survey area, ∆~Brandom. Note that pc stands forpadding cells.Method ∆Bxrandom ∆Byrandom ∆Bzrandom ∆Bxrandom ∆Byrandom ∆BzrandomFrequency of 1 HzMFV+Arithmetic 8.24 8.56 10.90MFV+Geometric 6.80 6.15 10.34MFV+Harmonic 16.48 16.06 22.33MSFV 3.30 3.34 4.61MSFV+O (2 pc) 0.90 0.86 0.26MSFV+O (4 pc) 1.13 1.04 0.62MSFV+O (8 pc) 1.26 1.18 0.86Frequency of 10 Hz Frequency of 100 HzMFV+Arithmetic 7.78 8.08 10.33 2.90 3.09 4.12MFV+Geometric 6.58 5.95 10.07 2.82 2.45 5.43MFV+Harmonic 16.08 15.66 21.88 7.38 7.00 11.40MSFV 3.14 3.18 4.41 1.22 1.25 1.87MSFV+O (2 pc) 0.95 0.91 0.38 0.94 0.89 0.57MSFV+O (4 pc) 1.17 1.09 0.71 1.06 0.99 0.76MSFV+O (8 pc) 1.30 1.22 0.92 1.06 0.98 0.82Frequency of 20 Hz Frequency of 200 HzMFV+Arithmetic 6.71 6.97 9.02 2.23 2.44 2.79MFV+Geometric 6.02 5.43 9.38 1.76 1.57 3.73MFV+Harmonic 15.03 14.62 20.69 4.81 4.60 7.69MSFV 2.76 2.78 3.92 0.97 1.04 1.26MSFV+O (2 pc) 1.04 1.00 0.56 0.78 0.75 0.45MSFV+O (4 pc) 1.26 1.18 0.86 0.84 0.79 0.58MSFV+O (8 pc) 1.36 1.28 1.04 0.83 0.77 0.61Frequency of 40 Hz Frequency of 400 HzMFV+Arithmetic 4.72 4.91 6.56 1.83 2.04 1.76MFV+Geometric 4.70 4.21 7.76 0.92 0.95 2.21MFV+Harmonic 12.23 11.83 17.42 2.68 2.68 4.18MSFV 1.98 1.99 2.92 0.83 0.92 0.71MSFV+O (2 pc) 1.12 1.07 0.71 0.62 0.63 0.35MSFV+O (4 pc) 1.31 1.24 0.99 0.64 0.63 0.43MSFV+O (8 pc) 1.37 1.30 1.12 0.62 0.60 0.44114Table 4.6: Average run time in seconds per simulation required to compute∆~Brandom on a two hexa-core Intel Xeon X5660 CPUs at 2.8 Hz with 64GB shared RAM using MATLAB. Setup time: time required to computethe local interpolation matrices and to assemble the reduced system ofequations to be solved. Solve time: time to solve the reduced system ofequations. Total time: the sum of the setup and solve times.Method Setup time (s) Solve time (s) Total time (s)MFV on fine OcTree mesh - 721 s 721 sMFV + Arithmetic - 114 s 114 sMFV + Geometric - 115 s 115 sMFV + Harmonic - 115 s 115 sMSFV 42 s 32 s 74 sMSFV+O (1 padding cells) 74 s 32 s 106 sMSFV+O (2 padding cells) 124 s 32 s 156 sMSFV+O (3 padding cells) 406 s 33 s 439 s1154.5 SummaryIn this chapter, we developed an oversampling technique for the multiscale FVmethod proposed in [73] to simulate quasi-static EM responses in the frequencydomain for geophysical settings that include highly heterogeneous conductive me-dia and features varying at different spatial scales. Simulating these types of geo-physical settings requires both large CPU time and memory usage; they oftenneed a very large and fine mesh to be discretized accurately. To the best of myknowledge, this is the first investigation focused on developing an oversamplingtechnique for geophysical EM problems in the multiscale literature.The method begins by assuming a coarse mesh nested into a fine mesh,which accurately discretizes the geophysical setting. For each coarse cell, weindependently solve a local version of the original Maxwell system subject to linearboundary conditions on an extended domain. To solve the local Maxwell system,we use the fine mesh contained in the extended domain and the MFV method(although a FE approach can be used for this part as well). Afterwards, these localsolutions, called basis functions, together with a weak-continuity condition areused within a Galerkin approach to construct a coarse-mesh (projected) versionof the global problem that is much cheaper to solve.For the examples presented, the proposed oversampling method significantlyimproves the accuracy relative to the MSFV method (without oversampling) pro-posed in [73] at the cost of more computing power, which is still lower than thecost of the fine-mesh solution. Our method produces results comparable to thoseobtained by simulating EM responses using MFV on a fine mesh, while drasticallyreducing the size of the linear system of equations and the computational time.Using the oversampling technique in the presented examples in combination withan OcTree mesh (i.e., an adaptive mesh refinement technique), the size of thecoarse-mesh system is only about 10% of the fine-mesh system size, while therelative error is significantly reduced for all the cases considered.Although the examples presented demonstrate a significant advantage of us-ing the MSFV with oversampling (as opposed to using MSFV without oversam-pling) to forward modeling geophysical EM responses in highly heterogeneousconductive media, there may be some cases where the proposed oversampling116technique may fail. For example, when the underlying conductivity structure con-tains features with large connectivities (e.g. a steel-cased well). For such cases,one alternative is to use large oversampling domains, which may increase signifi-cantly the cost of the forward simulation.117Chapter 5Concluding RemarksA smooth sea never made a skillful sailor. — Franklin D. Roosevelt5.1 Summary of Work and Research HighlightsThe aim of this dissertation is to investigate the applicability and feasibility of up-scaling and multiscale methods for efficiently modeling geophysical quasi-staticEM frequency-dependent problems. Modeling EM responses in complex geo-physical settings is crucial to the exploration, imaging and characterization ofburied natural resources, such as mineral and ground water deposits, and hy-drocarbon reserves. However, in practice, accurately simulating these types ofgeophysical settings can be computationally very expensive.Realistic EM geophysical settings typically consider large computational do-mains, features that vary at multiple spatial scales, and a wide variation over sev-eral orders of magnitude of the physical properties of the heterogeneous media(e.g. electrical conductivity). Since all of these factors can have a significant im-pact on the behavior of the EM responses of interest, if we need to obtain anaccurate approximation to the EM responses, the simulation mesh used in tra-ditional discretization techniques (e.g. FE or FV) should be able to capture thestructure of the heterogeneity present in the setting with sufficient detail. Thisneed leads to the use of a very large mesh to discretize the model, which resultsinto solving a very large, and often very ill-conditioned, system of equations —118in some cases, in the order of millions, or even larger than billions of unknowns.Such a large system of equations require specialized computing resources (e.g.clusters) to be solved. When an EM simulation is conducted in practice (e.g. fordifferent frequencies or survey configurations) or within an EM inversion proce-dure (where there is the need to deal with adjoint operators per frequency and/orsource), several forward EM simulations must be conducted [68]. This can lead toa very computational expensive process overall, and therefore it is of interest toreduce effectively the computational cost of individual EM simulations.In satisfying the aim of this investigation, I initially focused on identifying thecore mathematical ideas of some of the successful upscaling and multiscale meth-ods for the problem of modeling single-phase fluid flow in highly heterogeneousporous media, whose leading mathematical model is given by a linear elliptic PDE(i.e., Poisson’s equation), which is a real and scalar PDE. Such a problem sharesseveral key challenges similar to the problem of simulating geophysical EM re-sponses in highly heterogeneous settings, namely, the governing PDE model inboth problems is linear, the simulation considers large-scale computational do-mains, features varying at multiple spatial scales, and a wide variation over severalorders of magnitude of the physical properties (e.g. permeability) of the media.According to the literature review conducted (see Chapter 2), there is extensiveevidence demonstrating that for linear flow in porous media problems, upscalingand multiscale FE/FV methods are able to drastically reduce the cost of the sim-ulation, while yielding results whose accuracy is similar to that of those obtainedwith FE or FV on a fine mesh. Surprisingly, when I began this investigation, littleresearch was done along these lines for the problem of simulating geophysicalEM responses in highly heterogeneous settings.Once the core mathematical ideas of some of the successful upscaling andmultiscale FE/FV methods developed for flow in porous media problems wereidentified (see Section2.6), I extended these ideas to develop two new simulationmethods for the quasi-static Maxwell’s equations in the frequency domain:1. An upscaling framework for the electrical conductivity, which is thoroughlydiscussed in Chapter 3. The investigation for this topic resulted in an ex-panded abstract [29] and a peer-reviewed publication [32].1192. A multiscale FV with oversampling method, which is thoroughly discussedin Chapter 4. The investigation for this topic resulted in two expanded ab-stracts [30] and [31], and a peer-reviewed publication [33].Both methods were implemented using the parallel computing toolbox of MAT-LAB, and they can use both tensor and OcTree meshes. In order to test theperformance of each method, field-inspired and synthetic examples that includehighly heterogeneous conductive media were considered.The examples presented demonstrate that both proposed upscaling and mul-tiscale with oversampling methods can be feasible and applicable to geophysicalquasi-static EM problems in the frequency domain. In addition, this investigationshows that the combination of both methods with locally refined OcTree meshesis particularly advantageous to tackle complex geophysical EM modeling prob-lems. This combination enables both proposed methods to drastically reduce thesize of the problem when using a large domain and a mesh that must capture thespatial distribution of the media heterogeneity outside the region where the EMresponses are measured. In practice, these two modeling aspects are part of theprinciple reasons why the cost of geophysical EM simulations is computationallyexpensive.The examples presented indicate that the accuracy of a coarse-mesh solu-tion obtained with the two methods proposed in this dissertation depends on twofactors. First, the given fine-mesh (conductivity) model should be accurately dis-cretized. The accuracy of the coarse-mesh solution is as good as the accuracy ofthe conventional discretization method used (in this case MFV) to solve the fine-mesh problem. Second, as it was shown in Sections 3.5 and 4.4, the coarse-meshneeds to be carefully designed because it can significantly impact the quality of thesolution. For the examples presented, the mesh is designed to maintain the sameresolution as the fine mesh where the EM responses are expected to vary rapidlyinside the survey area where they are measured (i.e., at the source(s) and re-ceivers locations). The estimate of the proper cell sizes of the coarse mesh insidethe survey area is based on the skin depth, which is a typical mesh-design crite-rion used in practice for geophysical EM forward modeling simulations [68, 138].Furthermore, for the conductivity models used to test the two proposed methods,120it is possible to aggressively coarsen the mesh outside the survey area withoutlosing much accuracy. However, aggressively coarsening the mesh outside thesurvey area may not be possible for all cases. For example, when the underlyingelectrical conductivity structure have large connectivities.The set of local Dirichlet boundary conditions (shown in Table 3.3) used tosolve the local problems inside every coarse-mesh cell for both proposed upscal-ing and multiscale with oversampling approaches provided reasonable estimatesfor the examples presented in this dissertation. However, this set of boundaryconditions may not be appropriate to compute accurate coarse-mesh approxima-tions using the proposed two methods for all cases. In this investigation, the useof an extended local domain to solve the local problems mitigated the effect of im-posing the chosen local boundary conditions, which may cause a reduction on theaccuracy of the coarse-mesh solution for highly heterogeneous conductive media.For example, such effect can be particularly observed when high contrast of con-ductivity is present at the boundary of a coarse-mesh cell. Although using a localextended domain significantly improved the accuracy of the coarse-mesh solutionfor the cases considered, there may be other cases where this approach may notbe sufficient. For example, when the underlying conductivity structure containsfeatures with rather long connectivities (e.g. a steel-cased well). For such cases,one alternative is to use larger local extended domains, which may come at thecost of more computing time. Similar observations regarding the effect of localboundary conditions and the use of local extended domains are consistent withwhat has been reported in the upscaling and multiscale literature for flow in porousmedia problems.The proposed upscaling method requires a more cautious setup to performefficiently. It involves a number of parameters to be decided upon depending onthe purpose and requirements of the simulation. If the involved parameters arewell chosen, the proposed upscaling method retrieves accurate results. Other-wise, the proposed upscaling method can also be quite inaccurate. In general,upscaling methods lack the build-in feature to transition from the coarse mesh tothe fine mesh that multiscale approaches have and, often, they lack generality (cf.[23, 51]). The proposed upscaling method is not an exception to this observation.However, if there is the need to generate a coarse-scale electrical conductivity121model from a fine-scale conductivity model, as in the well log conductivity exam-ple presented in Section 3.3, the proposed upscaling method offers quite reliableresults to do so.The 3D examples presented in this dissertation, demonstrate that multiscalewith oversampling can retrieve accurate solutions at a fraction of the cost of thefine-mesh solution, even when the coarsening is extreme according to the empir-ical skin depth mesh-design criteria. In contrast, computing an upscaled solutionusing the proposed upscaling framework is significantly more expensive than com-puting a fine-mesh or a multiscaled solution. As it was shown in Chapter 3, thecomputational cost comes from solving a local parameter estimation problem forevery coarse cell. In addition, multiscale methods have the built-in feature to inter-polate the solution from the coarse mesh to the fine mesh, which opens the doorfor a natural extension of this methodology to a multilevel method. I will elaboratemore about such potential multilevel extension in Section 5.3.The next two sections summarize the work done and highlight the specificcontributions for each of the two methods proposed in this dissertation.5.1.1 Upscaling: Summary and HighlightsChapter 3 proposes an upscaling framework for the electrical conductivity thatposes upscaling as a least-squares parameter estimation problem. This upscal-ing framework constructs a coarse-mesh conductivity model by solving a non-linear optimization problem for each coarse-mesh cell. Such optimization prob-lem minimizes the difference between the user-chosen EM responses, which arecomputed using the fine- and coarse-scale quasi-static Maxwell’s equations. Theconstructed coarse-mesh conductivity models can be scalar or full SPD tensormodels, depending on the purpose and context of the simulation. Contrary toother least-squares upscaling formulations in the literature for flow in porous me-dia problems (Section 2.6.1), the formulation proposed in this investigation doesnot require regularization and it is practical to tackle 3D geophysical EM problemswith arbitrary electrical conductivity structures (for details in the discussion seeSection 3.4). For the 1D and 3D examples presented, the size of the coarse-mesh system solved was roughly 10% of the fine-mesh system size, while the122relative errors (in the secondary fields) were less than 5%. To the best of myknowledge, this investigation constitutes the first upscaling approach proposedfor geophysical EM problems in the literature.The proposed upscaling method is more expensive than MFV in a fine meshas (local) optimization problems need to be solved per coarse-mesh cell; however,since each local problem is formulated independently of the others, the cost canbe reduced by using a more efficient parallel implementation of the method on amore powerful machine.The most important lesson learned from this investigation is that different ex-periments require different upscaling criteria that result in different upscaled quan-tities. That is, for a given fine-scale conductivity structure, there is no uniqueupscaled model which completely describes it. The examples presented demon-strate that using the proposed upscaling method for different frequencies in thesurvey configuration leads to different upscaled models. Similarly changing thetype of EM responses to be matched in the upscaling criterion leads to differentupscaled conductivities.5.1.2 Multiscale with Oversampling: Summary and HighlightsChapter 4 proposes a MSFV with an oversampling method. Contrary to the caseof upscaling methods, where no previous research work along that line was donefor geophysical EM problems, the MSFV method for the quasi-static Maxwell’sequations was first proposed in [73]. In this investigation, I found that MSFVcould produce quite inaccurate solutions due to the effect of the local boundaryconditions used to construct the multiscale basis functions in each coarse-meshcell. In order to solve this issue, I extended the oversampling technique origi-nally proposed by [89], which has been quite effective to tackle a related issuefor single-phase flow in porous media problems, for application to geophysical EMproblems. In the proposed oversampling technique, the multiscale basis functionsare computing using an extended local domain. Doing so, the multiscale basisfunctions take into account structures outside the core coarse cell and the effectof the local boundary conditions is reduced. Then, the oversampled basis togetherwith a weak-continuity condition are used to construct a coarse-mesh version of123the global problem. To the best of my knowledge, this is the first investigationfocused on developing an oversampling technique for geophysical EM problemsin the literature.For the 3D examples presented, the proposed oversampling technique signif-icantly increases the accuracy of the MSFV method at the expense of more com-putational run time, which however, depending on the size of the extended domainchosen, can still be considerably lower compared to the time of the fine-mesh so-lution. Although the examples presented demonstrate a significant advantage ofusing the MSFV with oversampling (as opposed to using MSFV without oversam-pling) to forward modeling geophysical EM responses in highly heterogeneousconductive media, there may be some cases where the proposed oversamplingtechnique may fail. For example, when the underlying conductivity structure con-tains features with large connectivities (e.g. a steel-cased well). For such cases,one alternative is to use large oversampling domains, which may increase signifi-cantly the cost of the forward simulation.5.2 Impact to the Field of Computational Methods inGeophysical ElectromagneticsThe impact of this investigation to the current field of computational methods ingeophysical EM can be summarized in the following three outcomes.On a theoretical level, one outcome of this study is the advancement of knowl-edge by contributing two new methods to efficiently solve the quasi-static Max-well’s equations in the frequency domain with highly discontinuous coefficientsand features at multiple length scales: (1) an upscaling method and (2) a multi-scale FV with oversampling method. Both methods have the potential to be usefulfor other applications that have the same underlying mathematical model, for ex-ample, in the fields of medical imaging and the various engineering branches,such as electrical, mechanical or materials engineering.On a practical level, a second outcome of this study is that it lays the foun-dation to develop more efficient and usable parallel computing environments forboth multiscale and upscaling methods to simulate EM responses in complex geo-physical settings. In particular, these simulation environments combine upscaling124and multiscale methods with adaptive mesh refinement techniques. Doing so in-creases our current predictive and analytic capabilities by making the simulationof EM responses in larger and more complex geophysical settings more feasiblethan currently is possible.A third outcome of this study is that it demonstrates a novel use of upscalingand multiscale methods for geophysical EM applications. Developments in up-scaling and multiscale methods have been primarily within the field of modelingflow in porous media. This investigation extends the ideas developed for suchmethods for application to Maxwell’s equations, which is more involved due to thecomplex, vectorial nature of the Maxwell system.5.3 Remaining Challenges and Future ResearchThis study demonstrates the feasibility and applicability of both upscaling andmultiscale approaches to simulate EM responses in complex geophysical settings.Nonetheless, there remain some challenges to overcome and there are someinteresting paths to extend both research directions.A natural first extension for both upscaling and multiscale methods proposedhere is to adapt them for solving quasi-static EM problems in the time domain.In addition, both proposed methods can be adapted to use a different base dis-cretization method to solve the local Maxwell problems, such as a mimetic FD([121]) or FE ([104, 131], instead of the MFV method (Section 2.4).For both proposed upscaling and multiscale methods, the open question re-mains: ‘what is the best set of boundary conditions to solve the local problemson each coarse-mesh cell?’ An idea for selecting better local boundary conditionswould be to incorporate global information in the solution of the problem, as it isdone in the upscaling community for local-global approaches (e.g. [49, 61, 112]for flow in porous media upscaling methods.)This study assumes that only the electrical conductivity (Σ) captures the het-erogeneity of the geophysical setting. This is a reasonable assumption as in mostgeological environments, variations in the Earth’s magnetic permeability are in-significant and EM surveys are only sensitive to contrasts in the Earth’s electricalconductivity [138]. However, there are geophysical applications where the mag-125netic susceptibility (µ) can be heterogeneous as well. For example, when mod-eling settings that consider steel-well casing or some ore-bearing rocks, whichhave very high magnetic permeabilities. An extension of the upscaling frame-work for µ is possible using the least-squares formulations proposed in Chapter3. However, it will require a rigorous investigation to make appropriate choicesfor the number of parameters involved (e.g. data to be matched in the upscalingcriterion, boundary conditions for local problems, etc.). Contrary to the upscalingapproach, handling both heterogeneous µ and Σ simultaneously in the multiscaleFV with oversampling method proposed in Chapter 4 is relatively straight forward.The discretization of the local Maxwell systems required to compute the local pro-jection matrices in this case can be done as described in Section 2.4. The rest ofthe multiscale-related processes remain the same.Other possible extensions for the upscaling framework proposed in this studyinclude the construction of complex anisotropic upscaled conductivities, which isof relevance for geophysical EM applications where induced polarization effectsappear (e.g. [125]). In addition, the formulation can be adapted to choose the datato be matched in the upscaling criterion among the different electric or magneticfields or fluxes (i.e., ~E, ~H,~B or ~J), or some combination of them. I elaborated moreabout these ideas in Sections 3.2 and 3.4.Another interesting research path is to extend the two-level multiscale FVmethod for the quasi-static Maxwell’s equations into a multilevel method. One ofthe major challenges along this line is on how to construct a hierarchy of coarsespaces and their corresponding inter-level transfer operators, as well as the cor-responding approximation properties that guarantee proper bounds for the dis-cretization error. Once this first (major) challenge has been figure it out, the nextmajor step has to do with designing an efficient parallelization for such a mul-tilevel method. As discussed in Section 2.6.3, the multilevel multiscale mimeticmethod (M3) ([119]) and the element-based AMG (AMGe) method ([113, 140])have shown promising results along these lines for flow in porous media prob-lems. I have conducted some preliminary work along these two directions, whichI describe below.• This study lays the foundation for the construction of oversampled multi-126scale basis functions for the electric field, ~E (for every coarse-mesh cell)in a two-level setting. However, in order to construct a hierarchy of coarsespaces required for a robust mimetic multilevel extension, it is also of in-terest to generate the basis functions for ~B. In collaboration with somecolleagues at UBC, we are currently investigating the extension of the (two-level) M3 method to the Maxwell’s equations. This investigation focuses onconstructing mimetic multiscale basis functions for both ~E and ~B, while as-suming that both µ and Σ are highly discontinuous (heterogeneous). Theextension also considers both frequency- and time-domain EM problems.We have obtained some preliminary positive results. Such results are notdocumented in this dissertation, since they are out of the scope of this studyand they require a significant amount of effort to constitute a full furtherstudy.• While being a summer research intern at Lawrence Livermore National Lab-oratory, I conducted some preliminary studies on using their (C++) AMGelibrary to upscale the quasi-static Maxwell’s equations in the time domain.Such study focused on investigating the multilevel capability of the AMGe li-brary and its computational performance. In this case, I also obtained somepreliminary positive results. Such preliminary results are not documentedin this dissertation, since they are off-topic and they require a significantamount of effort to constitute a full further study. However, research alongthis direction seems promising to speed up the computational performancefor further geophysical EM problems.There are several other research directions to be explored that can lead tointeresting extensions for the two methods presented in this dissertation. Onlythe immediate extensions that I envision where discussed above; however, onecan continue to gain inspiration from the current successful developments andapplications of both upscaling and multiscale methods in the field of petroleumengineering.This investigation proposed two innovative mathematical alternatives to im-prove the computational performance of solving realistic geophysical EM forwardproblems in a more sustainable way than currently is possible. It also opened the127door to use such alternatives to create more powerful computational environmentscapable of simulating EM responses in larger and more complex geophysical set-tings than currently is possible. Working on advancing this type of computationaltools is important because they have a tangible impact in the overall processes ofimaging and monitoring buried natural resources using geophysical EM methods.128Bibliography[1] A. Abdulle, E. Weinan, B. Engquist, and E. Vanden-Eijnden. Theheterogeneous multiscale method. Acta Numer., 21(April 2012):1–87, apr2012. ISSN 0962-4929. doi:10.1017/S0962492912000025. → pages 34[2] R. Alcouffe, A. Brand, J. Dendy, and J. Painter. The Multi-grid method forthe difussion equation with strongly discontinuous coefficients. Siam J.Sci. Stat. Comput., 2(4):430–454, 1981.doi:http://epubs.siam.org/doi/abs/10.1137/0902035. → pages 85[3] P. Amestoy, I. Duff, J. Y. L’Excellent, and J. Koster. MUMPS: a generalpurpose distributed memory sparse solver. Int. Work. Appl. ParallelComput. Springer Berlin Heidelb., pages 121–130, 2001. → pages 29, 75,94[4] D. A. Aruliah and U. Ascher. Multigrid Preconditioning for Krylov Methodsfor Time-Harmonic Maxwell’s Equations in Three Dimensions. SIAM J. Sci.Comput., 24(2):702–718, jan 2002. ISSN 1064-8275.doi:10.1137/S1064827501387358. → pages 28[5] D. A. Aruliah, U. Ascher, E. Haber, and D. W. Oldenburg. A method for theforward modelling of 3-D electromagnetic quasi-static problems. Math.Model. Methods Appl. Sci., 11(1):1–21, 2001.doi:10.1142/S0218202501000702. → pages 13, 27, 28[6] U. Ascher. Numerical methods for evolutionary differential equations.Society for Industrial and Applied Mathematics, 2008. ISBN978-0-898716-52-8. → pages 8, 20, 24[7] D. Avdeev. Three-dimensional electromagnetic modelling and inversionfrom theory to application. Surv. Geophys., 26(6):767–799, 2005.doi:10.1007/s10712-005-1836-2. → pages 13129[8] D. Avdeev and S. Knizhnik. 3D integral equation modeling with a lineardependence on dimensions. Geophysics, 74(5):F89–F94, sep 2009. ISSN0016-8033. doi:10.1190/1.3190132. → pages 14[9] I. Babukat, G. Calozt, and J. E. Osborn. Special Finite Element Methodsfor a Class of Second Order Elliptic Problems with Rough Coefficients.SIAM J. Numer. ANAL, 31(4):945–981, 1994. ISSN 0036-1429.doi:10.1137/0731051. → pages 34, 35[10] I. Babusˇka and J. E. Osborn. Generalized finite element methods: Theirperformance and their relation to mixed methods. SIAM J. Numer. Anal.,20(3):510–536, 1983. ISSN 0036-1429. doi:10.1137/0720034. → pages35[11] I. Babuska, U. Banderjee, and J. E. Osborn. Survey of meshless andgeneralized finite element methods: A unified approach. Acta Numer., 12:1–125, 2003. ISSN 09624929. doi:10.1017/S0962492902000090. →pages 34[12] E. Badea, M. Everett, G. Newman, and O. Biro. Finite-element analysis ofcontrolled-source electromagnetic induction using Coulomb-gaugedpotentials. Geophysics, 66(3):786–799, 2001. ISSN 00168033.doi:10.1190/1.1444968. → pages 28[13] P. Benner and K. Willcox. A Survey of Projection-Based Model ReductionMethods for Parametric Dynamical Systems . SIAM Rev., 57(4):483–531,2015. → pages 35[14] A. Bensoussan, J.-L. Lions, and G. Papanicolau. Asymptotic Analysis forPeriodic Structures. North-Holland Pub. Co., New York, vol. 5 edition,1978. ISBN 0444851720. → pages 30[15] M. Benzi. Preconditioning Techniques for Large Linear Systems: ASurvey. J. Comput. Phys., 182(2):418–477, 2002. ISSN 00219991.doi:10.1006/jcph.2002.7176. → pages 28[16] A. J. Bergqvist and S. G. Engdahl. A homogenization procedure of fieldquantities in laminated electric steel. In IEEE Trans. Magn., volume 37,pages 3329–3331, 2001. doi:10.1109/20.952606. → pages 43[17] A. Bermu´dez, D. Go´mez, and P. Salgado. Eddy-current losses inlaminated cores and the computation of an equivalent conductivity. IEEETrans. Magn., 44(12):4730–4738, 2008. ISSN 00189464.doi:10.1109/TMAG.2008.2005118. → pages 43130[18] J. G. Berryman and G. M. Hoversten. Modelling electrical conductivity forearth media with macroscopic fluid-filled fractures. Geophys. Prospect., 61(2):471–493, mar 2013. ISSN 00168025.doi:10.1111/j.1365-2478.2012.01135.x. → pages 31[19] R. U. Bo¨rner. Numerical modelling in geo-electromagnetics: advances andchallenges. Surv. Geophys., 31(2):225–245, oct 2010. ISSN 0169-3298.doi:10.1007/s10712-009-9087-x. → pages 3, 4, 13[20] A. Borzı` and V. Schulz. Computational Optimization of Systems Governedby Partial Differential Equations - Chapter 1. Society for Industrial andApplied Mathematics, 2012. ISBN 9781611972054.doi:10.1137/1.9781611972054. → pages 46[21] A. Brandt. Multi-level adaptive solutions to boundary-value problems.Math. Comput., 31(138):333–390, 1977. ISSN 0025-5718.doi:10.1090/S0025-5718-1977-0431719-X. → pages 36[22] A. Brandt. Multiscale scientific computation: Review 2001. Multiscale andmultiresolution methods, (6682):3–95, 2002.doi:10.1007/978-3-642-56205-1. → pages 4, 6, 15, 36[23] A. Brandt. Principles of systematic upscaling. In Multiscale Methods Bridg.Scales Sci. Eng., chapter 7, pages 1–19. Oxford University Press, 2009.ISBN 9780191715532. doi:10.1093/acprof:oso/9780199233854.003.0007.→ pages 15, 121[24] A. Brandt and E. O. Livne. Multigrid Techniques: 1984 Guide withApplications to Fluid Dynamics, Revised Edition. Society for Industrial andApplied Mathematics, 2011. ISBN 978-1-61197-074-6.doi:http://dx.doi.org.ezproxy.library.ubc.ca/10.1137/1.9781611970753. →pages 37[25] R. P. Brent. Algorithms for Minimization Without Derivatives, volume 19.Prentice-Hall, 1972. ISBN 0130223352. doi:10.1109/TAC.1974.1100629.→ pages 50[26] M. Brezina, A. J. Cleary, R. D. Falgout, V. E. Henson, J. E. Jones, T. A.Manteuffel, S. F. McCormick, and J. W. Ruge. Algebraic Multigrid Basedon Element Interpolation (AMGe). SIAM J. Sci. Comput., 22(5):1570–1592, 2001. ISSN 1064-8275. doi:10.1137/S1064827598344303.→ pages 37, 38131[27] W. L. Briggs, V. E. Henson, and S. F. McCormick. A Multigrid Tutorial.Society for Industrial and Applied Mathematics, 2nd edition, 2000. ISBN9780898714623.doi:http://dx.doi.org.ezproxy.library.ubc.ca/10.1137/1.9780898719505. →pages 28[28] L. Cao, Y. Zhang, W. Allegretto, and Y. Lin. Multiscale Asymptotic Methodfor Maxwell’s Equations in Composite Materials. SIAM J. Numer. Anal., 47(6):4257–4289, 2010. doi:10.1137/080741276. → pages 31[29] L. A. Caudillo-Mata, E. Haber, L. J. Heagy, and D. W. Oldenburg.Numerical upscaling of electrical conductivity: A problem specificapproach to generate coarse-scale models, chapter 130, pages 680 – 684.Society of Exploration Geophysicists, 2014.doi:10.1190/segam2014-1488.1. → pages 40, 119[30] L. A. Caudillo-Mata, E. Haber, and C. Schwarzbach. An OversamplingTechnique for Multiscale Finite Volume Method to SimulateFrequency-domain Electromagnetic Responses. In SEG Tech. Progr.Expand. Abstr. 2016, pages 981 – 985, 2016.doi:10.1190/segam2016-13974625.1. → pages 82, 120[31] L. A. Caudillo-Mata, E. Haber, and C. Schwarzbach. A Multiscale FiniteVolume Method with Oversampling for Geophysical ElectromagneticSimulations. In EAGE - ECMOR XV - 15th Eur. Conf. Math. Oil Recover. -Expand. Abstr., 2016. doi:10.3997/2214-4609.201601892. → pages 82,120[32] L. A. Caudillo-Mata, E. Haber, L. J. Heagy, and C. Schwarzbach. Aframework for the upscaling of the electrical conductivity in the quasi-staticMaxwell’s equations. J. Comput. Appl. Math., 317:388–402, 2017. ISSN03770427. doi:10.1016/j.cam.2016.11.051. → pages 40, 119[33] L. A. Caudillo-Mata, E. Haber, and C. Schwarzbach. An oversamplingtechnique for the multiscale finite volume method to simulateelectromagnetic responses in the frequency domain. Comput. Geosci.,pages 1–18, 2017. ISSN 15731499. doi:10.1007/s10596-017-9647-y. →pages 82, 120[34] J. Chen and Q. H. Liu. Discontinuous Galerkin Time-Domain Methods forMultiscale Electromagnetic Simulations: A Review. Proc. IEEE, 101(2):242–254, feb 2013. ISSN 0018-9219. doi:10.1109/JPROC.2012.2219031.→ pages 13132[35] T. Chen. New Methods for Accurate Upscaling with Full-tensor Effects.PhD thesis, Stanford University, 2009. → pages 33, 34, 58, 64[36] T. Chen, M. G. Gerritsen, J. V. Lambers, and L. J. Durlofsky. Globalvariable compact multipoint methods for accurate upscaling with full-tensoreffects. Comput. Geosci., 14(1):65–81, 2010.doi:10.1007/s10596-009-9133-2. → pages 33, 64[37] Y. Chen, L. J. Durlofsky, M. G. Gerritsen, and X. H. Wen. A coupledlocalglobal upscaling approach for simulating flow in highly heterogeneousformations. Adv. Water Resour., 26(10):1041–1060, 2003.doi:10.1016/S0309-1708(03)00101-5. → pages 33[38] C. Christensen. Multilevel techniques for Reservoir Simulation. Doctoraldissertation, Technical University of Denmark, 2016. → pages 8, 38[39] E. Chung, Y. Efendiev, and T. Y. Hou. Adaptive multiscale model reductionwith Generalized Multiscale Finite Element Methods. J. Comput. Phys.,320:69–95, 2016. ISSN 10902716. doi:10.1016/j.jcp.2016.04.054. →pages 8, 36[40] E. T. Chung and J. Zou. The Eigenvalues and Eigenspaces of SomeDiscrete Div- and Curl-Related Operators. SIAM J. Matrix Anal. Appl., 24(4):1149, 2003. ISSN 08954798. doi:10.1137/S0895479801382483. →pages 27[41] M. Commer and G. A. Newman. An accelerated time domain finitedifference simulation scheme for three-dimensional transientelectromagnetic modeling using geometric multigrid concepts. Radio Sci.,41(3):1–15, 2006. ISSN 00486604. doi:10.1029/2005RS003413. →pages 14[42] M. Commer, M. B. Kowalsky, J. Doetsch, G. A. Newman, and S. Finsterle.MPiTOUGH2: A parallel parameter estimation framework for hydrologicaland hydrogeophysical applications. Comput. Geosci., 65:127–135, 2014.ISSN 00983004. doi:10.1016/j.cageo.2013.06.011. → pages 1[43] T. A. Davis. Direct Methods for Sparse Linear Systems. Society forIndustrial and Applied Mathematics, 1st edition, 2006. ISBN 0898716136,9780898716139. → pages 7, 29[44] J. Dendy. Black box multigrid. J. Comput. Phys., 48(3):366–386, dec 1982.ISSN 00219991. doi:10.1016/0021-9991(82)90057-2. → pages 85133[45] A. J. Desbarats. Numerical estimation of effective permeability insand-shale formations. Water Resour. Res., 23(2):273–286, 1987. ISSN19447973. doi:10.1029/WR023i002p00273. → pages 32[46] A. Druinsky, P. Ghysels, X. S. Li, O. Marques, S. Williams, A. Barker,D. Kalchev, and P. Vassilevski. Comparative Performance Analysis ofCoarse Solvers for Algebraic Multigrid on Multicore and ManycoreArchitectures, pages 116–127. Springer International Publishing, Cham,2016. ISBN 978-3-319-32149-3. doi:10.1007/978-3-319-32149-3 12. →pages 101[47] L. J. Durlofsky. Numerical calculation of equivalent grid block permeabilitytensors for heterogeneous porous media. Water Resour. Res., 27(5):699–708, may 1991. ISSN 00431397. doi:10.1029/91WR00107. → pages61[48] L. J. Durlofsky. Coarse scale models of two phase flow in heterogeneousreservoirs: volume averaged equations and their relationship to existingupscaling techniques. Comput. Geosci., 2(2):73–92, mar 1998. ISSN1573-1499. doi:10.1023/A:1011593901771. → pages 33, 56, 59[49] L. J. Durlofsky. Upscaling of geocellular models for reservoir flowsimulation: a review of recent progress. 7th Int. Forum Reserv. Simul.Bu¨hl/Baden-Baden, Ger., pages 1–58, 2003. → pages 8, 30, 31, 32, 33,34, 46, 51, 56, 58, 59, 61, 62, 64, 79, 125[50] L. J. Durlofsky, Y. Efendiev, and V. Ginting. An adaptive local-globalmultiscale finite volume element method for two-phase flow simulations.Adv. Water Resour., 30(3):576–588, 2007. ISSN 03091708.doi:10.1016/j.advwatres.2006.04.002. → pages 33, 36[51] Y. Efendiev and T. Y. Hou. Multiscale finite element methods: theory andapplications. Springer New York, 2009. ISBN 9780387094953.doi:10.1007/978-0-387-09496-0. → pages 8, 30, 34, 35, 36, 59, 85, 86,87, 88, 121[52] Y. Efendiev, J. Galvis, G. Li, and M. Presho. Generalized Multiscale FiniteElement Methods. Oversampling Strategies. Int. J. Multiscale Comput.Eng., 12(6):465–484, 2014. doi:10.1615/IntJMultCompEng.2014007646.→ pages 36[53] R. Eso and D. W. Oldenburg. 3D forward modeling and inversion of CSEMdata at the San Nicolas massive sulphide deposit, 2007. → pages 13134[54] R. Eymard, T. Galloue¨t, and R. Herbin. Finite volume methods (Schemesand Analysis course at the University of Wroclaw), 2000. → pages 14[55] R. D. Falgout, P. S. Vassilevski, and L. T. Zikatanov. On two-gridconvergence estimates. Numer. Linear Algebr. with Appl., 12(5-6):471–494, 2005. ISSN 10705325. doi:10.1002/nla.437. → pages 37[56] C. L. Farmer. Upscaling: a review. Int. J. Numer. Methods Fluids, 40(1-2):63–78, 2002. doi:10.1002/fld.267. → pages 8, 30, 31, 32, 34, 44, 46, 51,59, 61, 62, 64, 79[57] C. G. Farquharson and D. W. Oldenburg. Inversion of time-domainelectromagnetic data for a horizontally layered Earth. Geophys. J. Int., 114(3):433–442, 1993. ISSN 1365246X.doi:10.1111/j.1365-246X.1993.tb06977.x. → pages 50[58] D. Fleisch. A Student’s Guide to Maxwell’s Equations. CambridgeUniversity Press, 2008. ISBN 9780521877619. → pages 4, 13, 17, 19[59] K. Gao, S. Fu, R. L. Gibson, E. T. Chung, and Y. Efendiev. GeneralizedMultiscale Finite-Element Method (GMsFEM) for elastic wave propagationin heterogeneous, anisotropic media. J. Comput. Phys., 295:161–188,2015. ISSN 10902716. doi:10.1016/j.jcp.2015.03.068. → pages 8[60] M. G. Gerritsen and L. J. Durlofsky. Modeling fluid flow in oil reservoirs.Annu. Rev. Fluid Mech., 37(1):211–238, jan 2005. ISSN 0066-4189.doi:10.1146/annurev.fluid.37.061903.175748. → pages 8, 30, 31, 32, 33,61[61] M. G. Gerritsen and J. V. Lambers. Integration of localglobal upscaling andgrid adaptivity for simulation of subsurface flow in heterogeneousformations. Comput. Geosci., 12(2):193–208, 2008.doi:10.1007/s10596-007-9078-2. → pages 33, 125[62] P. Ghysels, X. S. Li, F. H. Rouet, S. Williams, and A. Napov. An efficientmulti-core implementation of a novel HSS-structured multifrontal solverusing randomized sampling. SIAM J. Sci. Comput., 38(5):1–26, 2015.ISSN 10957200. doi:10.1137/15M1010117. → pages 29[63] G. H. Golub and C. Greif. On Solving Block-Structured Indefinite LinearSystems. SIAM J. Sci. Comput., 24(6):2076–2092, 2003. ISSN1064-8275. doi:10.1137/S1064827500375096. → pages 7, 28135[64] C. Greif and D. Scho¨tzau. Preconditioners for the discretizedtime-harmonic maxwell equations in mixed form. Numer. Linear Algebr.with Appl., 14(4):281–297, 2007. → pages 7, 28[65] A. Gupta and H. Avron. WSMP: Watson Sparse Matrix Package, Part I-direct solution of symmetric systems. Technical Report 98462, IBMCorporation, 2017. → pages 29[66] E. Haber. A mixed finite element method for the solution of themagnetostatic problem with highly discontinuous coefficients in 3D.Comput. Geosci., 4(4):323–336, 2000. doi:10.1023/A:1011540222718. →pages 28[67] E. Haber. A parallel method for large scale time domain electromagneticinverse problems. Appl. Numer. Math., 58(4):422–434, apr 2008. ISSN01689274. doi:10.1016/j.apnum.2007.01.017. → pages 46[68] E. Haber. Computational Methods in Geophysical Electromagnetics.Society for Industrial and Applied Mathematics, 2014. ISBN978-1-611973-79-2. → pages 3, 4, 5, 6, 7, 13, 14, 15, 16, 17, 19, 21, 22,24, 25, 26, 27, 28, 29, 46, 49, 62, 64, 67, 73, 76, 119, 120[69] E. Haber. Advanced Modeling and Inversion Tools for Controlled SourceEM. In 76th EAGE Conf. Exhib., number June 2014, pages WS9–B02.EAGE, 2014. doi:10.3997/2214-4609.20140555. → pages 85, 86[70] E. Haber and U. Ascher. Fast finite volume simulation of 3Delectromagnetic problems with highly discontinuous coefficients. SIAM J.Sci. Comput., 22(6):1943–1961, jan 2001. ISSN 1064-8275.doi:10.1137/S1064827599360741. → pages 7, 13, 22, 28[71] E. Haber and U. Ascher. Preconditioned all-at-once methods for large ,sparse parameter estimation problems. Inverse Probl., 17(6):1847–1864,2001. doi:10.1088/0266-5611/17/6/319. → pages 7[72] E. Haber and S. Heldmann. An octree multigrid method for quasi-staticMaxwell’s equations with highly discontinuous coefficients. J. Comput.Phys., 223(2):783–796, may 2007. ISSN 00219991.doi:10.1016/j.jcp.2006.10.012. → pages 7, 26, 75, 76, 94, 105[73] E. Haber and L. Ruthotto. A multiscale finite volume method for Maxwell’sequations at low frequencies. Geophys. J. Int., 199(2):1268 – 1277, 2014.doi:10.1093/gji/ggu268. → pages 11, 36, 81, 82, 83, 85, 86, 96, 101, 116,123136[74] E. Haber and L. Ruthotto. A multiscale finite volume method for Maxwell’sequations at low frequencies. Geophys. J. Int., 199(2):1268–1277, 2014.doi:10.1093/gji/ggu268. → pages 34[75] E. Haber, U. Ascher, D. A. Aruliah, and D. W. Oldenburg. Fast simulation of3D electromagnetic problems using potentials. J. Comput. Phys., 163(1):150–171, 2000. doi:10.1006/jcph.2000.6545. → pages 27[76] E. Haber, U. Ascher, and D. W. Oldenburg. 3D forward modelling of timedomain electromagnetic data. 2002 SEG Annu. Meet., (1):1–4, 2002. →pages 13[77] R. Haberman. Elementary Applied Partial Differential Equations WithFourier Series and Boundary Value Problems, 1987. → pages 13[78] H. Haddar, R. Hiptmair, P. Monk, and R. Rodriguez. ComputationalElectromagnetism. Springer International Publishing, Cetraro, Italy, 2015edition, 2015. ISBN 9783319193052. doi:10.1007/978-3-319-19306-9. →pages 4[79] H. Hajibeygi and P. Jenny. Multiscale finite-volume method for parabolicproblems arising from compressible multiphase flow in porous media. J.Comput. Phys., 228:5129–5147, 2009. ISSN 00219991.doi:10.1016/j.jcp.2009.04.017. → pages 35, 87[80] H. Hajibeygi and P. Jenny. Adaptive iterative multiscale finite volumemethod. J. Comput. Phys., 230(3):628–643, feb 2011. ISSN 00219991.doi:10.1016/j.jcp.2010.10.009. → pages 8[81] H. Hajibeygi, G. Bonfigli, M. Hesse, and P. Jenny. Iterative multiscalefinite-volume method. J. Comput. Phys., 227(19):8604–8621, oct 2008.ISSN 00219991. doi:10.1016/j.jcp.2008.06.013. → pages 8, 35[82] P. Henning and D. Peterseim. Oversampling for the multiscale finiteelement method. SIAM Multiscale Model. Simul., 11(4):1149–1175, 2013.doi:10.1137/120900332. → pages 36[83] R. Hiptmair. Multigrid Method for Maxwell’s Equations. SIAM J. Numer.Anal., 36(1):204–225, 1998. doi:10.1137/S0036142997326203. → pages7, 28[84] R. Hiptmair. Analysis of Multilevel Methods for Eddy Current Problems.AMS Math. Comput., 72(243):1281–1303, 2003.doi:10.1090/S0025-5718-02-01468-0. → pages 28137[85] L. Holden and B. F. Nielsen. Global upscaling of permeability inheterogeneous reservoirs; the Output Least Squares (OLS) method.Transp. Porous Media, 40(2):115–143, 2000. ISSN 01693913.doi:10.1023/A:1006657515753. → pages 34, 46, 61[86] L. Horesh and E. Haber. A Second Order Discretization of Maxwell’sEquations in the Quasi-Static Regime on OcTree Grids. SIAM J. Sci.Comput., 33(5):2805–2819, 2011. doi:10.1137/100798508. → pages 7,26, 75, 76, 94, 105[87] M. F. Horstemeyer. Practical Aspects of Computational Chemistry.Springer Netherlands, Dordrecht, 2010. ISBN 978-90-481-2686-6.doi:10.1007/978-90-481-2687-3. → pages 8[88] T. Y. Hou. Numerical Approximations to Multiscale Solutions in PartialDifferential Equations. In Front. Numer. Anal., pages 241–301. Springer,2003. → pages 8, 87[89] T. Y. Hou and X. Wu. A Multiscale Finite Element Method for EllipticProblems in Composite Materials and Porous Media. J. Comput. Phys.,134(1):169–189, jun 1997. ISSN 00219991. doi:10.1006/jcph.1997.5682.→ pages 8, 34, 35, 36, 83, 85, 87, 88, 90, 123[90] T. Y. Hou, X. Wu, and Z. Cai. Convergence of a Multiscale Finite ElementMethod for Elliptic Problems with Rapidly Oscillating Coefficients. Math.Comput. Am. Math. Soc., 68(227):913–943, 1999.doi:10.1155/2013/764165. → pages 35, 36[91] L. J. Hughes. A Short History of Electrical Techniques in PetroleumExploration. Technical report, Zonge Engineering & ResearchOrganization, Inc., Tucson, Arizona, U.S.A., 2013. → pages 1[92] T. J. Hughes, G. R. Feijo´o, L. Mazzei, and J. B. Quincy. The variationalmultiscale methoda paradigm for computational mechanics. Comput.Methods Appl. Mech. Eng., 166(1-2):3–24, 1998. ISSN 00457825.doi:10.1016/S0045-7825(98)00079-6. → pages 34[93] J. M. Hyman and M. Shashkov. Adjoint operators for the naturaldiscretizations of the divergence, gradient and curl on logically rectangulargrids. Appl. Numer. Math., 25(4):413–442, 1997.doi:10.1016/S0168-9274(97)00097-4. → pages 15, 17, 18, 19, 21138[94] J. M. Hyman and M. Shashkov. Natural Discretizations for the Divergence,Gradient, and Curl on Logically Rectangular Grids. Comput. Math. withAppl., 33(4):81–104, 1997. doi:10.1016/S0898-1221(97)00009-6. →pages 18, 19, 21[95] J. M. Hyman and M. Shashkov. Mimetic discretizations for Maxwellequations and the equations of magnetic diffusion. Technical report,United States. Department of Energy. Office of Energy Research, 1998. →pages 21[96] J. M. Hyman and M. Shashkov. Mimetic discretizations for Maxwell’sequations. J. Comput. Phys., 151(2):881–909, may 1999. ISSN 00219991.doi:10.1006/jcph.1999.6225. → pages 15, 21, 25[97] J. M. Hyman and M. Shashkov. The orthogonal decomposition theoremsfor mimetic finite difference methods. SIAM J. Numer. Anal., 36(3):788–818, jan 1999. ISSN 0036-1429. doi:10.1137/S0036142996314044.→ pages 21[98] J. M. Hyman and M. Shashkov. Mimetic Finite Difference Methods forMaxwell’s Equations and the Equations of Magnetic Diffusion. Prog.Electromagn. Res., 32:89–121, 2001. → pages 17, 18[99] U. S. Inan and R. A. Marshall. Numerical Electromagnetics: The FDTDMethod. Cambridge University Press, 2011. ISBN 9780521190695. →pages 14[100] H. Jahandari, S. Ansari, and C. G. Farquharson. Comparison betweenstaggered grid finitevolume and edgebased finiteelement modelling ofgeophysical electromagnetic data on unstructured grids. J. Appl.Geophys., 138:185–197, 2017. doi:10.1016/j.jappgeo.2017.01.016. →pages 6, 14, 28[101] P. Jenny, S. Lee, and H. Tchelepi. Multi-scale finite-volume method forelliptic problems in subsurface flow simulation. J. Comput. Phys., 187(1):47–67, may 2003. ISSN 00219991. doi:10.1016/S0021-9991(03)00075-5.→ pages 8, 34, 35, 36, 83, 87[102] X. Jiang and W. Zheng. Homogenization of Quasi-Static Maxwell’sequations. Multiscale Model. Simul., 12(1):152–180, 2014.doi:10.1137/130921398. → pages 31139[103] V. V. Jikov, O. A. Oleinik, and S. M. Kozlov. Homogenization of differentialoperators and integral functionals. Springer-Verlag, 1st edition, 1994.ISBN 9783642846618. doi:10.1007/978-3-642-84659-5. → pages 30, 31[104] J. Jin. The Finite Element Method in Electromagnetics. Wiley, New York,2nd edition, 2002. ISBN 978-0-471-43818-2. → pages 14, 15, 17, 18, 21,22, 59, 84, 85, 125[105] J. E. Jones and B. Lee. A multigrid method for variable coefficientmaxwell’s equations. SIAM J. Sci. Comput., 27(5):1689–1708, 2006.doi:10.1137/040608283. → pages 28[106] A. Kelbert, N. Meqbel, G. D. Egbert, and K. Tandon. ModEM: A modularsystem for inversion of electromagnetic geophysical data. Comput.Geosci., 66:40–53, 2014. ISSN 00983004.doi:10.1016/j.cageo.2014.01.010. → pages 13, 14[107] C. T. Kelley. Iterative methods for optimization. Society for Industrial andApplied Mathematics, 1999. ISBN 0898714338, 9780898714333.doi:http://dx.doi.org.ezproxy.library.ubc.ca/10.1137/1.9781611970920. →pages 62[108] K. Key and J. Ovall. A parallel goal-oriented adaptive finite elementmethod for 2.5D electromagnetic modelling. Geophys. J. Int., 186(1):137–154, 2011. doi:10.1111/j.1365-246X.2011.05025.x. → pages 7[109] S. Knapek. Matrix-dependent multigrid homogenization for diffusionproblems. SIAM J. Sci. Comput., 20(2):515—-533 (electronic), 1998. ISSN1064-8275. doi:10.1137/S1064827596304848. → pages 37[110] G. Kristensson and N. Wellander. Homogenization of the Maxwellequations at fixed frequency. SIAM J. Appl. Math., 64(1):170–195, 2003.ISSN 0036-1399. doi:10.1137/S0036139902403366. → pages 31, 43[111] A. Kumar, L. Farmer, G. L. Jerauld, and D. Li. Efficient upscaling fromcores to simualtion models. In SPE Annu. Tech. Conf. Exhib.SPE-38744-MS, pages 1–16, San Antonio, Texas, 1997. Society ofPetroleum Engineers. doi:10.2118/38744-MS. → pages 30[112] J. V. Lambers, M. G. Gerritsen, and D. Fragola. Multiphase 3-D FlowSimulation with Integrated Upscaling, MPFA Discretization and Adaptivity.In SPE Reserv. Simul. Symp. 2-4 February, Woodlands, Texas, numberSPE-118983-MS. Society of Petroleum Engineers, 2009.doi:10.2118/118983-MS. → pages 33, 125140[113] I. V. Lashuk and P. S. Vassilevski. The construction of the coarse de Rhamcomplexes with improved approximation properties. Comput. MethodsAppl. Math., 14(2):257–303, 2014. ISSN 16099389.doi:10.1515/cmam-2014-0004. → pages 38, 126[114] C. S. Lee. Generalization of Mixed Multiscale Finite Element Methods withApplications. Doctoral dissertation, Texas A&M University, 2016. → pages8, 36, 38[115] J. Li, C. G. Farquharson, and X. Hu. 3D vector finite-elementelectromagnetic forward modeling for large loop sources using a total-fieldalgorithm and unstructured tetrahedral grids. Geophysics, 82(1):E1–E16,2017. ISSN 0016-8033. doi:10.1190/geo2016-0004.1. → pages 14, 28[116] X. S. Li and J. W. Demmel. SuperLU DIST: A scalable distributed-memorysparse direct solver for unsymmetric linear systems. ACM Trans. Math.Softw., 29(2):110–140, 2003. ISSN 00983500.doi:10.1145/779359.779361. → pages 29[117] C. J. Lin and J. J. More´. Newton’s method for large bound-constrainedoptimization problems. SIAM J. Optim., 9(4):1100–1127, jan 1999. ISSN1052-6234. doi:10.1137/S1052623498345075. → pages 62[118] K. Lipnikov, J. Morel, and M. Shashkov. Mimetic finite difference methodsfor diffusion equations on non-orthogonal non-conformal meshes. J.Comput. Phys., 199:589–597, 2004. ISSN 00219991.doi:10.1016/j.jcp.2004.02.016. → pages 7[119] K. Lipnikov, J. D. Moulton, and D. Svyatskiy. A multilevel multiscale mimetic(M3) method for two-phase flows in porous media. J. Comput. Phys., 227:6727–6753, 2008. ISSN 00219991. doi:10.1016/j.jcp.2008.03.029. →pages 30, 38, 126[120] K. Lipnikov, J. D. Moulton, and D. Svyatskiy. Adaptive Strategies in theMultilevel Multiscale Mimetic Method for Two-Phase Flows in PorousMedia. Multiscale Model. Simul., 9(3):991–1016, 2011. ISSN 1540-3459.doi:10.1137/100787544. → pages 38[121] K. Lipnikov, G. Manzini, and M. Shashkov. Mimetic finite differencemethod. J. Comput. Phys., 257(PB):1163–1227, 2014. ISSN 10902716.doi:10.1016/j.jcp.2013.07.031. → pages 15, 18, 125[122] S. MacLachlan. Improving Robustness in Multiscale Methods. PhD thesis,University of Colorado, 2004. → pages 8, 15, 30, 31, 32, 39141[123] S. MacLachlan and J. D. Moulton. Multilevel upscaling through variationalcoarsening. Water Resour. Res., 42(2):W02418:1–9, 2006.doi:10.1029/2005WR003940.1. → pages 8, 36, 39, 83[124] N. K. Madsen and R. W. Ziolkowski. A Three-Dimensional Modified FiniteVolume Technique for Maxwell’s Equations. Electromagnetics, 10(1-2):147–161, 1990. doi:10.1080/02726349008908233. → pages 25[125] D. Marchant, E. Haber, L. Beran, and D. W. Oldenburg. 3D modeling of IPeffects on electromagnetic data in the time domain. SEG 2012 Expand.Abstr., pages 1–5, 2012. → pages 126[126] A. Menshov, Y. Brick, C. Torres-Verdin, and A. E. Yilmaz. Recent progressin rigorous algorithms for the fast solution of 3-D EM frequency-domainintegral-equations. In 6th Int. Symp. Three-Dimensional Electromagn.Berkeley,, pages 1–4, Berkely, California, 2017. → pages 14[127] G. W. Milton. The theory of composites. Cambridge University Press,Cambridge; New York, 2002. ISBN 0521781256. → pages 31[128] Y. Mitsuhata and T. Uchida. 3D magnetotelluric modeling using the T-Ωfinite-element method, 2004. ISSN 00168033. → pages 28[129] R. Mittra. Computational Electromagnetics Recent Advances andEngineering Applications. Springer New York, 2014. ISBN9781461443810. → pages 4[130] P. Monk. An analysis of Nedelec’s method for the spatial discretization ofMaxwell’s equations. J. Comput. Appl. Math., 47(1):101–121, 1993.doi:10.1016/0377-0427(93)90093-Q. → pages 14, 90[131] P. Monk. Finite element methods for Maxwell’s equations. Clarendon,2003. ISBN 9780198508885, 0198508883.doi:10.1093/acprof:oso/9780198508885.001.0001. → pages 14, 15, 18,19, 21, 22, 25, 59, 84, 85, 90, 125[132] J. D. Moulton, J. E. Dendy, and J. M. Hyman. The black box multigridnumerical homogenization algorithm. J. Comput. Phys., 108(142):80–108,1998. doi:10.1006/jcph.1998.5911. → pages 8, 32[133] W. A. Mulder. Geophysical modelling of 3D electromagnetic diffusion withmultigrid. Comput. Vis. Sci., 11(3):129–138, 2008.doi:10.1007/s00791-007-0064-y. → pages 28142[134] G. A. Newman. A Review of High-Performance Computational Strategiesfor Modeling and Imaging of Electromagnetic Induction Data, 2014. ISSN01693298. → pages 3, 4, 6, 7, 13, 27, 28[135] G. A. Newman and D. L. Alumbaugh. Frequency-domain modelling ofairborne electromagnetic responses using staggered finite differences.Geophys. Prospect., 43(8):1021–1042, nov 1995. ISSN 0016-8025.doi:10.1111/j.1365-2478.1995.tb00294.x. → pages 18, 28[136] B. F. Nielsen and A. Tveito. An upscaling method for one-phase flow inheterogeneous reservoirs . A weighted output least squares ( WOLS )approach. Comput. Geosci., 2:93–123, 1998. ISSN 14200597.doi:10.1023/A:1011541917701. → pages 34, 44, 46, 61[137] J. Nocedal and S. J. Wright. Numerical Optimization. Springer New York,2nd edition, 2006. ISBN 9780387303031.doi:10.1007/978-0-387-40065-5. → pages 62, 67[138] D. W. Oldenburg and Y. Li. Inversion for applied geophysics: a tutorial. InNear-Surface Geophys., chapter 5, pages 89–150. Society of ExplorationGeophysicists, 2005. ISBN 978-1-56080-130-6.doi:10.1190/1.9781560801719.ch5. → pages 3, 9, 120, 125[139] D. W. Oldenburg and D. A. Pratt. Geophysical inversion for mineralexploration: A decade of progress in theory and practice. Proc. Explor.,2007. → pages 1, 3, 4[140] J. E. Pasciak and P. S. Vassilevski. Exact de Rham sequences on spacesdefined on macro-elements in two and three spatial dimensions. SIAM J.Sci. Comput., 30(5):2427–2446, 2008.doi:https://doi.org/10.1137/070698178. → pages 38, 126[141] G. Pavliotis and A. Stuart. Multiscale methods: averaging andhomogenization. Springer New York, 2008. ISBN 9780387738284.doi:10.1007/978-0-387-73829-1. → pages 8, 31[142] A. Raiche. Modelling and inversion progress, problems and challenges.Surv. Geophys., 15(2):159–207, 1994. doi:10.1007/BF00689859. →pages 13[143] P. Renard and G. de Marsily. Calculating equivalent permeability: a review.Adv. Water Resour., 20(5-6):253–278, 1997. ISSN 03091708.doi:10.1016/S0309-1708(96)00050-4. → pages 30143[144] J. M. Reynolds. An introduction to applied and environmental geophysics.John Wiley, 1997. ISBN 9780471955559, 0471955558, 9780471968023,0471968021. → pages 1, 2, 7, 9, 16, 29, 48, 49[145] K. F. Riley, M. P. Hobson, and S. J. Bence. Mathematical methods forphysics and engineering. Cambridge University Press, 3rd edition, 2006.ISBN 0521861535, 9780521861533, 0521679710, 9780521679718. →pages 20, 21, 28[146] Y. Saad. Iterative Methods for Sparse Linear Systems. Society forIndustrial and Applied Mathematics, 2nd edition, 2003. ISBN9780898718003. doi:10.2277/0898715342. → pages 7, 27, 28[147] X. Sa´nchez-Vila, J. P. Girardi, and J. Carrera. A synthesis of approachesto upscaling of hydraulic conductivities. Water Resour. Res., 31(4):867–882, 1995. doi:10.1029/94WR02754. → pages 30[148] O. Schenk, K. Ga¨rtner, W. Fichtner, and A. Stricker. PARDISO: ahigh-performance serial and parallel sparse linear solver in semi-conductordevice simulation. Futur. Gener. Comput. Syst., 18(1):69–78, 2001. →pages 29[149] C. Schwarzbach. Stability of finite element discretization of Maxwell’sequations for geophysical applications. PhD thesis, University of Frieberg,Germany, 2009. → pages 7, 14, 15, 18, 28, 29[150] C. Schwarzbach and E. Haber. Finite element based inversion fortime-harmonic electromagnetic problems. Geophys. J. Int., 193(2):615–634, feb 2013. ISSN 0956-540X. doi:10.1093/gji/ggt006. → pages 14[151] B. Shafiro and M. Kachanov. Anisotropic effective conductivity of materialswith nonrandomly oriented inclusions of diverse ellipsoidal shapes. J.Appl. Phys., 87(12):8561–8569, 2000. ISSN 00218979.doi:10.1063/1.373579. → pages 31[152] J. Smith. Conservative modeling of 3-D electromagnetic fields, Part I:Properties and error analysis. Geophysics, 61(5):1308–1318, sep 1996.ISSN 0016-8033. doi:10.1190/1.1444054. → pages 28[153] R. Smith. Electromagnetic Induction Methods in Mining Geophysics from2008 to 2012. Surv. Geophys., 35(1):123–156, apr 2013. ISSN0169-3298. doi:10.1007/s10712-013-9227-1. → pages 1, 3, 4, 13144[154] L. P. Song, S. D. Billings, L. R. Pasion, and D. W. Oldenburg. TransientElectromagnetic Scattering of a Metallic Object Buried in UnderwaterSediments. IEEE Trans. Geosci. Remote Sens., 54(2):1091–1102, 2016.ISSN 01962892. doi:10.1109/TGRS.2015.2473851. → pages 1[155] K. Steklova and E. Haber. Joint hydrogeophysical inversion: stateestimation for seawater intrusion models in 3D. Comput. Geosci., 21(1):75–94, 2017. ISSN 15731499. doi:10.1007/s10596-016-9595-y. → pages1[156] J. A. Stratton. Electromagnetic Theory. John Wiley & Sons, Inc., 2007. →pages 4, 13, 16[157] R. Streich. 3D finite-difference frequency-domain of controlled-sourceelectromagnetic data: Direct solution and optimization for high accuracy.Geophysics, 74(5):F95–F105, 2009. doi:10.1190/1.3196241. → pages 13,28[158] A. Tarantola. Logarithmic parameters. 2001. → pages 13[159] A. Tarantola. Inverse problem theory and methods for model parameterestimation. Society for Industrial and Applied Mathematics, Paris, 2005.ISBN 0898715725, 9780898715729. → pages 3, 4[160] H. Tchelepi, P. Jenny, S. Lee, and C. Wolfsteiner. An Adaptive MultiphaseMultiscale Finite Volume Simulator for Heterogeneous Reservoirs. In SPEReserv. Simul. Symp., pages 1–10, 2005. doi:10.2118/93395-MS. →pages 35[161] W. M. Telford, L. P. Geldart, and R. E. Sheriff. Applied Geophysics, volume2nd. Cambridge University Press, 2nd edition, 1990. ISBN9781139167932. doi:http://dx.doi.org/10.1017/CBO9781139167932.009.→ pages 1, 2, 5, 9, 16[162] The MathWorks Inc. MATLAB and Statsistics Toolbox Release 2011b.Technical report, Massachusetts, United States, 2014. → pages 50, 75, 94[163] S. Torquato. Random heterogeneous materials: microstructure andmacroscopic properties. Springer, 2002. ISBN 0387951679. → pages 30,31, 67, 109[164] U. Trottenberg, C. Oosterlee, and A. Schuller. Multigrid. Academic Press,2000. ISBN 012701070X. → pages 28, 37145[165] E. S. Um, M. Commer, G. A. Newman, and G. M. Hoversten. Finiteelement modelling of transient electromagnetic fields near steel-casedwells. Geophys. J. Int., 202(2):901–913, 2015. ISSN 0956-540X.doi:10.1093/gji/ggv193. → pages 18[166] M. A. Valle´e, R. S. Smith, and P. Keating. Metalliferous mining geophysicsState of the art after a decade in the new millennium. Geophysics, 76(4):W31–W50, 2011. doi:10.1190/1.3587224. → pages 1, 3, 4, 13[167] P. S. Vassile. Coarse Spaces By Algebraic Multigrid: MultigridConvergence and Upscaling Error Estimates. Adv. Adapt. Data Anal., 03(01n02):229–249, 2011. ISSN 1793-5369.doi:10.1142/S1793536911000830. → pages 36, 37, 38[168] P. E. Wannamaker and G. W. Hohmann. Electromagnetic modeling ofthree-dimensional bodies in layered earths using integral equations.Technical Report January, University of Utah Research Institute, EarthScience Laboratory, 1982. → pages 14, 28[169] S. H. Ward and G. W. Hohmann. Electromagnetic theory for geophysicalapplications. In M. N. Nabighian, editor, Electromagn. methods Appl.Geophys., volume 1, chapter 4, pages 130–311. Society of ExplorationGeophysicists, 1988. ISBN 978-1-56080-263-1. → pages 4, 9, 13, 14, 16,17, 18, 26, 42, 46, 49, 73, 88[170] J. E. Warren and H. S. Price. Flow in Heterogeneous Porous Media. Soc.Pet. Eng. J., 1(03):153–169, 1961. ISSN 0197-7520. doi:10.2118/1579-G.→ pages 32[171] E. Weinan, B. Engquist, X. Li, W. Ren, and E. Vanden-Eijnden.Heterogeneous multiscale methods: a review. Commun. Comput. Phys., 2(3):367–450, 2007. → pages 34[172] C. J. Weiss and G. A. Newman. Electromagnetic induction in a generalized3D anisotropic earth, Part 2: The LIN preconditioner. Geophysics, 68(3):922–930, 2003. ISSN 00168033. doi:10.1190/1.1581044. → pages 7, 28[173] N. Wellander. Homogenization of the Maxwell’s equations: case I. LinearTheory. Appl. Math., 46(1):29–51, 2001. doi:10.1023/A:1013727504393.→ pages 31, 43[174] N. Wellander. Homogenization of the Maxwell equations: Case II.Nonlinear conductivity. Appl. Math., 47(3):255–283, 2002.doi:10.1023/A:1021797505024. → pages 31, 43146[175] X. H. Wen and J. J. Go´mez-Herna´ndez. Upscaling hydraulic conductivitiesin heterogeneous media: An overview. J. Hydrol., 183(1):ix – xxxii, 1996.doi:10.1016/S0022-1694(96)80030-8. → pages 8, 30, 31, 32, 33, 34, 59,61, 79[176] D. Wynne, M. Attalla, T. Berezniuk, H. Berhane, D. K. Cotterill, R. Strobl,and D. M. Wightman. Athabasca Oil Sands Database.McMurray/Wabiskaw Deposit. Alberta Geological Survey. Technical report,Alberta Research Council, Edmonton, 1994. → pages 49[177] D. Yang, D. Fournier, and D. W. Oldenburg. 3D Inversion of EM data atLalor mine: in pursuit of a unified electrical conductivity model. In Explor.Deep VMS Ore Bodies Hudbay Lalor Case Study, pages 1–4. BCGeophysical Society, 2014. → pages 71, 92[178] K. Yee. Numerical Solution of Initial Boundary Value Problems InvolvingMaxwell’s Equations in Isotropic Media. Antennas Propagation, IEEETrans., 14(3):302–307, 1966. → pages 14, 18[179] M. S. Zhdanov. Electromagnetic geophysics: Notes from the past and theroad ahead. Geophysics, 75(5):75A49–75A66, sep 2010. ISSN0016-8033. doi:10.1190/1.3483901. → pages 1, 3, 4, 9, 13147
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Multiscale and upscaling methods for geophysical electromagnetic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Multiscale and upscaling methods for geophysical electromagnetic forward modeling Caudillo Mata, Luz Angélica 2017
pdf
Page Metadata
Item Metadata
Title | Multiscale and upscaling methods for geophysical electromagnetic forward modeling |
Creator |
Caudillo Mata, Luz Angélica |
Publisher | University of British Columbia |
Date Issued | 2017 |
Description | Accurate and efficient simulation of electromagnetic responses in realistic geophysical settings is crucial to the exploration, imaging, and characterization of buried natural resources, such as mineral and hydrocarbon deposits. However, in practice, these simulations are computationally expensive. The geophysical settings consider highly heterogeneous media and features at multiple spatial scales that require a very large mesh to be accurately represented. This results in a system of equations to be solved that often exceeds the limits of average computers. Thus, the key is to reduce the problem size but retain the accuracy of the electromagnetic responses. Upscaling and multiscale techniques have been successfully applied to the problem of simulating fluid flow through heterogeneous porous media, where they are able to drastically reduce the size of the resulting fine-mesh system by casting it into a coarse-mesh system that is much cheaper to solve, while achieving a level of accuracy similar to that obtained with conventional discretization schemes. Recognizing the success that such techniques have had in fluid flow applications, this dissertation extends their use for application to electromagnetic modeling. In this dissertation, two new parallel simulation methods for the quasi-static Maxwell’s equations in the frequency domain are proposed: an upscaling framework for the electrical conductivity, and a multiscale finite volume with oversampling method. Both methods are combined with an adaptive mesh refinement technique (OcTree) to boost their computational performance. The performance of these methods is demonstrated by using field-inspired and synthetic examples that include a large electrical conductivity contrast. This investigation shows that both proposed methods are feasible to tackle geophysical electromagnetic problems, where being able to reduce the size of the problem can be particularly advantageous when extended domains are considered or when the mesh must capture the spatial distribution of the media heterogeneity outside the region where the electromagnetic responses are measured. Furthermore, both methods are new contributions to the literature in the field of computational methods in geophysical electromagnetics. Finally, both methods increase the current predictive and analytic capabilities by making the simulation of electromagnetic responses in larger and more complex geophysical settings more feasible than currently is possible. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2017-10-18 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0357172 |
URI | http://hdl.handle.net/2429/63343 |
Degree |
Doctor of Philosophy - PhD |
Program |
Geophysics |
Affiliation |
Science, Faculty of Earth, Ocean and Atmospheric Sciences, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2017-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2017_november_caudillomata_luzangelica.pdf [ 8.08MB ]
- Metadata
- JSON: 24-1.0357172.json
- JSON-LD: 24-1.0357172-ld.json
- RDF/XML (Pretty): 24-1.0357172-rdf.xml
- RDF/JSON: 24-1.0357172-rdf.json
- Turtle: 24-1.0357172-turtle.txt
- N-Triples: 24-1.0357172-rdf-ntriples.txt
- Original Record: 24-1.0357172-source.json
- Full Text
- 24-1.0357172-fulltext.txt
- Citation
- 24-1.0357172.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.24.1-0357172/manifest