# Eldad Haber

#### Relevant Degree Programs

## Graduate Student Supervision

##### Doctoral Student Supervision (Jan 2008 - May 2019)

Inverse modeling is a powerful tool for extracting information about the subsurface from geophysical and hydrologic data. Geophysical inverse problems are inherently multidisciplinary, requiring elements from the relevant physics, numerical simulation, and optimization, as well as knowledge of the geologic setting, hydrologic processes, and a comprehension of the interplay between all of these elements. Increasingly geoscientists are tackling complex problems that require integration of multiple types of information in order to better characterize the subsurface. However, many of the sub-fields of geophysics are developing simulation and inversion approaches, algorithms, and supporting software in isolation. This isolation is a barrier to quantitative integration and leads to inefficiencies in advancing interdisciplinary research. Greater efficiencies, and higher quality outcomes, could be achieved if (hydro)geophysicists had a common framework to accelerate an integrated approach. The main goal of my thesis is to organize the components of (hydro)geophysical simulations and inverse problems, and synthesize these into a comprehensive, modular framework.The development of a geophysical framework requires considering a number of disciplines and geophysical problems (e.g. electromagnetics and potential fields) to ensure generality as well as extensibility. However, the goal is also to have the framework work outside of geophysics and most notably in hydrogeology; vadose zone fluid flow is used as a model problem. Fluid flow in the vadose zone is governed by the Richards equation; it is parameterized by hydraulic conductivity, which is a nonlinear function of pressure head. The computational scalability of the Richards equation inversion is a significant challenge for three dimensional inversions in hydrogeophysics. Existing work explicitly calculates the sensitivity matrix using finite difference or automatic differentiation, however, for large-scale problems these methods are constrained by computation and memory. This dissertation provides an implicit sensitivity algorithm that enables large-scale inversion problems for distributed parameters in the Richards equation to become tractable on modest computational resources.

View record

In the modern era of diminishing returns on fixed exploration budgets, challenging targets, and ever-increasing numbers of multi-parameter datasets, proper management and integration of available data is a crucial component of any mineral exploration program. Machine learning algorithms have successfully been used for years by the technology sector to accomplish just this task on their databases, and recent developments aim at appropriating these successes to the field of mineral exploration. Framing the exploration task as a supervised learning problem, the geological, geochemical and geophysical information can be used as training data, and known mineral occurences can be used as training labels. The goal is to parameterize the complex relationships between the data and the labels such that mineral potential can be estimated in under-explored regions using available geoscience data. Numerous models and algorithms have been attempted for mineral prospectivity mapping in the past, and in this thesis we propose two new approaches. The first is a modified support vector machine algorithm which incorporates uncertainties on both the data and the labels. Due to the nature of geoscience data and the characteristics of the mineral prospectivity mapping problem, uncertainties are known to be very important. The algorithm is demonstrated on a synthetic dataset to highlight this importance, and then used to generate a prospectivity map for copper-gold porphyry targets in central British Columbia using the QUEST dataset as a case study. The second approach, convolutional neural networks, was selected due to its inherent sensitivity to spatial patterns. Though neural networks have been used for mineral prospectivity mapping, convolutional neural nets have yet to be applied to the problem. Having gained extreme popularity in the computer vision field for tasks involving image segmentation, identification and anomaly detection, the algorithm is ideally suited to handle the mineral prospectivity mapping problem. A CNN code is developed in Julia, then tested on a synthetic example to illustrate its effectiveness at identifying coincident structures in a multi-modal dataset. Finally, a subset of the QUEST dataset is used to generate a prospectivity map using CNNs.

View record

Parameter and state estimation for groundwater models within a coupled hydrogeophysical framework has become common in the last few years as it has been shown that such estimates are usually better than those from a single data inversion. Different approaches have been suggested in literature to combine the essentially two different modalities in order to obtain better estimates for groundwater models, and improve monitoring of processes such as solute transport. However, the coupled approaches usually come at a price of higher computational cost and difficulties in coupling the geophysical and groundwater inverse problems. Unlike in other studies, we developed both the groundwater and geophysical models in the same computational environment in order to test different minimization strategies. When solving the coupled inverse problem, the objective function consists of data misfit and regularization terms as well as a coupling term that relates groundwater and geophysical states. We present a novel approach to solve the inverse problem using an Alternating Direction Method of Multipliers (ADMM) to minimize the coupled objective function. ADMM enables us to treat the groundwater and geophysical part separately and thus use existing software with minor changes.However, ADMM as well as many other coupled approaches relies on implementing some petrophysical relationship to couple the groundwater and geophysical variable. Such relationships are usually uncertain and hard to parametrize for a large region and can potentially produce solute mass errors in the final model estimates. Therefore, in this thesis we examine coupled approaches that replace the fixed petrophysical relationship by a more loose structure similarity constraint. Besides, we propose efficient computational methods to minimize the objective function when there is no explicit petrophysical constraint. All approaches were tested on 3D synthetic examples. In the solute tracer test we estimated hydraulic conductivity or solute distribution using a structure coupled inversion, and were able to reduce the errors compared to a single data inversion alone. For a more complex example of seawater intrusion we implemented the ADMM method, and obtained better estimates for the solute distribution compared to just considering each data separately, or solving the problem with a simple coupled approach.

View record

Airborne electromagnetic (AEM) data is commonly collected for detecting buried natural resources, and this technique is sensitive to subsurface electrical resistivity distributions. The subsequent process of 3D AEM inversion constructs a physical property model from this data in order to better understand the size and shape of potential hidden resources. This thesis is designed to develop practical strategies to improve 3D AEM inversion accuracy in geologic settings where AEM data sets produce inconsistent or unsatisfactory inversion results.In this research, two overarching problematic scenarios are examined. First, in regions where an AEM survey overlaps with other electromagnetic data sets, a novel cooperative approach is introduced. This method is first tested on synthetic data where instead of producing an inversion model from each data set, the cooperative algorithm finds one consistent 3D resistivity model with improved resolution. The approach is then applied to field data over a high-sulfidation epithermal gold deposit where similar improvements are observed.The second scenario relates to improving 3D AEM inversions over thin conductive anomalies, a common geophysical target for copper and gold deposits. A new parametric inversion is developed using skewed Gaussian ellipsoids to represent target bodies. The approach is general but applied to frequency and time-domain AEM data with one or multiple anomalies. Combined with a voxel algorithm, the parametric inversion forms a hybrid approach capable of resolving thin conductive targets with smooth surrounding features. This hybrid technique is tested on synthetic data over conductive targets in a resistive background, and consistently produces models that are easier to interpret compared to voxel inversions alone. Field examples from a volcanogenic massive sulfide and an orogenic gold deposit highlight the practical nature of this method to image conductive mineralization with increased precision.The thesis concludes by analyzing a setting where multiple spatially overlapping AEM data sets exist over narrow conductive anomalies. Here, parametric, cooperative and voxel inversions are applied together to generate a consistent 3D resistivity model with thin targets and smooth background features. This section includes a discussion about potential pitfalls of such an approach when incompatible field measurements are encountered.

View record

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.

View record

This thesis examines the derivation and implementation of the discrete adjoint method for several time-stepping methods. Our results are important for gradient-based numerical optimization in the context of large-scale model calibration problems that are constrained by nonlinear time-dependent PDEs. To this end, we discuss finding the gradient and the action of the Hessian of the data misfit function with respect to three sets of parameters: model parameters, source parameters and the initial condition. We also discuss the closely related topic of computing the action of the sensitivity matrix on a vector, which is required when performing a sensitivity analysis. The gradient and Hessian of the data misfit function with respect to these parameters requires the derivatives of the misfit with respect to the simulated data, and we give the procedures for computing these derivatives for several data misfit functions that are of use in seismic imaging and elsewhere. The methods we consider can be divided into two categories, linear multistep (LM) methods and Runge-Kutta (RK) methods, and several variants of these are discussed. Regular LM and RK methods can be used for ODE systems arising from the semi-discretization of general nonlinear time-dependent PDEs, whereas implicit-explicit and staggered variants can be applied when the PDE has a more specialized form. Exponential time-differencing RK methods are also discussed. The implementation of the associated adjoint time-stepping methods is discussed in detail. Our motivation is the application of the discrete adjoint method to high-order time-stepping methods, but the approach taken here does not exclude lower-order methods. All of the algorithms have been implemented in MATLAB using an object-oriented design and are written with extensibility in mind. For exponential RK methods it is illustrated numerically that the adjoint methods have the same order of accuracy as their corresponding forward methods, and for linear PDEs we give a simple proof that this must always be the case. The applicability of some of the methods developed here to pattern formation problems is demonstrated using the Swift-Hohenberg model.

View record

Imaging and prediction of fluid flow within the subsurface provides information crucial to decision making processes in fields such as groundwater management and enhanced oil recovery. The flow of a fluid through a reservoir depends primarily on the permeability of the subsurface rock; a quantity that is often unknown throughout the entire domain of the reservoir. One method for predicting flow is to estimate the permeability of the reservoir and simulate flow through a mathematical subsurface flow model. Given the model, flow data can be inverted to estimate the permeability. However, this inversion approach can lead to inaccurate results due to the sparse sampling of flow data, and thus inaccurate predictions. To acquire a higher sampling of data, geophysical survey techniques are applied in order to efficiently collect a higher density of data sampled at the surface. These data are sensitive to changes to the geophysical properties of the reservoir due to flow. Inversion of geophysical data then provides images of changes to the geophysical properties of the reservoir. In order to estimate the flow parameters using geophysical data, the two mathematical models require coupling. The thesis therefore proposes two approaches to improve the imaging and prediction of flow. First, a novel coupled inverse problem for estimating the fluid velocity field and the initial geophysical property model from geophysical data is developed. Second, a new method of optimally designing the geophysical survey for the coupled inverse problem is developed. The new adaptive design approach builds on traditional A-Optimal design methods such that historic data are included in the design algorithm. This produces designs that adapt with flow in the subsurface and reduce the collection of unnecessary data. Both the coupled inverse problem and adaptive survey design method are demonstrated using a seismic tomography geophysical survey and a tracer advection fluid flow model. Numerical examples show that the coupled approach yields an improved flow estimate as well as improved image quality, while the adaptive optimal designs provide sufficient geophysical data.

View record

In this thesis, we introduce a new class of embarrassingly parallel parameter learning algorithms for Markov random fields (MRFs) with untied parameters, which are efficient for a large class of practical models. The algorithms parallelize naturally over cliques and, for graphs of bounded degree, have complexity that is linear in the number of cliques. We refer to these algorithms with the acronym LAP, which stands for Linear And Parallel. Unlike their competitors, the marginal versions of the proposed algorithms are fully parallel and for log-linear models they are also data efficient, requiring only the local sufficient statistics of the data to estimate parameters. LAP algorithms are ideal for parameter learning in big graphs and big data applications.The correctness of the newly proposed algorithms relies heavily on the existence and uniqueness of the normalized potential representation of an MRF. We capitalize on this theoretical result to develop a new theory of correctness and consistency of LAP estimators corresponding to different local graph neighbourhoods. This theory also establishes a general condition on composite likelihood decompositions of MRFs that guarantees the global consistency of distributed estimators, provided the local estimators are consistent. We introduce a conditional variant of LAP that enables us to attack parameter estimation of fully-observed models of arbitrary connectivity, including fully connected Boltzmann distributions. Once again, we show consistency for this distributed estimator, and relate it to distributed pseudo-likelihood estimators.Finally, for linear and non-linear inverse problems with a sparse forward operator, we present a new algorithm, named iLAP, which decomposes the inverse problem into a set of smaller dimensional inverse problems that can be solved independently. This parallel estimation strategy is also memory efficient.

View record

##### Master's Student Supervision (2010 - 2018)

Optimization problems with PDE constraints arise in many applications of science and engineering. The PDEs often describe physical systems and are characterized by a set of parameters, such as the material properties of the system, or velocity of the wave traveling through the earth’s subsurface. Usually the parameters cannot be directly measured and can only be inferred by solving inverse problems. Solving large scale inverse problems usually require prohibitively large number of PDE solves. As a result, various dimensionality reduction methods have been proposed to reduce the number of PDE solves. This work builds on a type of dimensionality method called the Simultaneous Source (SS) method. This method relies on random sampling and the stochastic trace estimator to reduce the number of PDE solves to a more manageable size. One of the limitations of the SS method is it can only be applied to data sets with no missing entries. Unfortunately, data sets with missing entries are common in practice, such as in geophysical applications.In this thesis, we propose a new coupled optimization problem that extend the SS method to data with missing entries. The block coordinate descent method (BCDM) is used to break the coupled problem into two subproblems. The first subproblem optimizes over the data and fills the missing entries by the rank minimization technique. The second subproblem minimizes over the model parameter and uses SS method to reduce the number of PDE solves into a more manageable size. We test our method in the context of a full-waveform inversion (FWI) problem and present our numerical results in which we are able to reduce the number of PDE solves by 98% with less than 2% difference between the model using 1 PDE solve and the model using full PDE solves.

View record

In this work, we implement an inversion algorithm for airborne electromagnetic (AEM) data in the frequency domain by using 2D conductivity models. First, we discretize the 2.5D Maxwell's equations on a staggered grid and test the numerical accuracy of the forward solution. The inverse problem is then solved by regularized minimization approach using the limited memory BFGS variant of the quasi-Newton method. Next, EM responses from a synthetic 2D conductivity model are inverted to validate the algorithm. Finally, we use the algorithm on an AEM field dataset from a RESOLVE survey and compare the inversion results to those obtained from a well-established 1D implementation.

View record

In this work, we implement an inversion algorithm for airborne electromagnetic (AEM) data in the frequency domain by using 2D conductivity models. First, we discretize the 2.5D Maxwell's equations on a staggered grid and test the numerical accuracy of the forward solution. The inverse problem is then solved by regularized minimization approach using the limited memory BFGS variant of the quasi-Newton method. Next, EM responses from a synthetic 2D conductivity model are inverted to validate the algorithm. Finally, we use the algorithm on an AEM field dataset from a RESOLVE survey and compare the inversion results to those obtained from a well-established 1D implementation.

View record

If this is your researcher profile you can log in to the Faculty & Staff portal to update your details and provide recruitment preferences.