Tutorial: Geophysical Inversion in SimPEG

all right hi everyone it's wonderful to have you all here I'm extremely excited to be participating in the transform virtual conference this year and I'm thrilled to be sharing with you some of the work that our team has been doing on building up sim pick a software package for simulation and parameter estimation in geophysics and so today we'll be walking through a few examples and I'm grateful you're all here so you'll see at the bottom of the screen there there's two links so this first link here this is where you can go to access the github repository that contains the notebooks and materials for this course and then there's also a link for the slack channel and it's this this hash tag here this channel here SOT 20 Tuesday simple is where you can go to ask questions and there's a whole team of folks there who are willing to willing to help answer so if you run into any issues please feel free to post there I so before we get going I wanted to say a big thank you to the organizers this has been a really exciting event to be a part of I think it's amazing that we can connect with folks from all around the world and be sharing some interesting work on the geosciences but in particular I want to thank Matt and Brendan who have been coordinating with me and created this opportunity to be presenting um so again here's the here's the link for their resources and this will this will be posted in a few different places and it's also available in the slack channel perfect so let's start by considering what what type of experiments are we going to be looking at and what is the set up of our problem that we're considering so here what we're going to be looking at is geophysical simulations and inversions and so we've got the earth which is characterized by some physical properties so maybe electrical resistivity and what we're going to do is we're going to input some source energy and measure some data that are sensitive to those variations and physical properties and then the goal in an inversion is to take those data and then try and obtain some information about the subsurface physical properties and so there's a whole bunch of different physical properties that we consider we can consider things like density magnetic susceptibility electrical properties like resistivity or charge ability or permittivity as well as elastic properties so in this tutorial we're going to focus in on electrical resistivity and charge ability and but a lot of the fundamental principles that were we exploring are going to be the same for all of these all of these different properties so let's start by considering electrical resistivity or conductivity and so what is that it's a measure of how easily current can pass through a material so there's two different symbols we often use to characterize this so Sigma we refer to as electrical conductivity which is in units of Siemens per meter and then there's resistivity which is its inverse so 1 over Sigma and that's in units of ohm meters and so we'll end up using these somewhat interchangeably I will do my best to stick with resistivity is the primary property that will be we'll be speaking about today but if I do say conductivity just know that it's the inverse of resistivity and so electrical resistivity it depends on quite a number of factors and we're looking at earth materials it's things like mineralogy the porosity and permeability and maybe the nature of the pore fluid so if we have a saline fluid versus a more fresh fluid that's gonna change the resistivity so here on the right we've just got a sample of a few different materials that we might encounter in geologic settings and you can see that a resistivity it varies over several orders of magnitude so we can go from anything like a massive sulfide which is quite conductive a down to igneous rocks which are very resistive and anything in between the other physical property who will consider is called chargeability and so what this is is it's the tendency for material to act like a capacitor and build up charges and so the actual physical mechanisms of this are quite they're involved in there's active research going on and really understanding the fundamental principles of what's going on but conceptually there's two models that are quite constructive to use a sort of a mental understanding of what's going on so the first that will consider is called membrane polarization and so what this is is if we have a poor throat that's narrow and narrow enough that we're not going to get ion ions in the fluid going through there if we apply an electric field we'll get a buildup of charges because they can't simply pass they can't simply pass through there and so that gives us a polarization and that's what we're sensitive to when we talk about a material being chargeable the other mechanism is called electrode polarization so if instead of having just a narrow poor throat we actually have some sort of metallic mineral in that poor throat then what happens when we apply an electric field we get we polarized that we polarized that metallic particle and then that's gonna attract charges and ions in the fluid on either side and so again that gives us this polarization and when we talk about earth materials so some of the materials that tend to be more chargeable are materials like pyrite or sulfides and so those are things that we are we're quite sensitive to when we do induced polarization surveys which are sensitive to charge ability and we'll talk about those in a moment okay so what's the basic geophysical experiment that we consider when we're looking for conductive or chargeable materials so what we're going to do is we're going to consider a mineral exploration example and so this is based on a deposit called the oil era deposit where we've got a conductive region of mineralization and it's capped by a resistive resistive unit here so what we're going to do is we have current electrodes that were going to connect to the ground and what we're going to do is we're going to inject current through those and then our receivers they're going to be these potential electrodes and they're going to measure a potential difference so first things first we turn on our electrodes and we set up the current so the currents in the earth they're gonna preferentially flow through conductive units so if here our mineralized units are more conductive than the surrounding house rock we're going to get a preferential flow of currents through those units and with done we then get charges that build up at interfaces so when we transition from a conductor to a resistor or a resistor to a conductor we get a buildup of charges at those interfaces and then as a results those charges create electric potentials and that is then what we measure at the surface of the earth so by measuring potentials at the surface of the earth we can get some information about where charges are building up and not is connected to then the resistivity of the subsurface so the first part that I walk through is typical of what we would call ADC resistivity experiment a direct current resistivity experiment so we turn on the curtains sorry we turn on our transmitter and those create currents that preferentially flow through conductors or around resistors and then we get the charge buildup and in a DC experiment that happens all instantaneously is that we get this buildup of charges instantaneously when we talk about induced polarization that charge buildup occurs over time and so in one experiment we can actually get information about both of these physical properties so here when we turn on a transmitter we instantaneously get a potential difference and then this slow buildup is the charging up or polarizing it materials and then when we shut that off we again the DC component goes away right away and then this decay is related to the the charge of the charge ability of those materials and so by looking at both the instantaneous response as well as these the case we can then get information both about the resistivity of the subsurface and the charge ability of the subsurface so the example that we're gonna look at today is in mineral exploration so DC and IP have been recognized as potentially useful for mineral exploration since the 1950s but it wasn't really until the 90s that we started to have the algorithms as well as the computational power to actually start working with these data and in order to recover 2d images of the subsurface and so why this is important is we actually want to try and understand and build up geologic images and geologic interpretations from our geophysical data so the specific deposit that we're gonna look at is called the century deposit so it is near the Mount eisah region in Queensland Australia and so these are the data that we're gonna be considering it is a zinc led silver deposit and this was really one of the early successes of DC resistivity and IP in mineral in mineral exploration as we'll see in the examples coming up IEP was really influential and helping understand the extent of mineralization in in this region and here there's a paper that we'll be referencing extensively this mutton 2000 paper this is linked here and I believe it should also be linked in the main resources that we shared on the github repository okay so let's start now by looking at what what's the geology and what are the physical properties that we're going to be looking looking for so the mineralized sequence that we're interested in it's about 40 meters thick and it's within these black carbonaceous shales so if we look here at this diagram if these four units they're tracking basically right right along here in these carbonaceous shales and so the two physical property is what we're going to be getting from both of these is the resistivity gives us structural information so here we see we've got some sand stones and siltstones and shales and so those might be more conductive units where we have a higher porosity for example and we can see a fault here where we've got a discontinuity between where these units where these units are in the subsurface and so that's gonna be the first step that we're gonna work with is working with these resistivity data and then the church ability is actually associated with the mineralized units and so in this case there's sulfides and so that's something that we know we can be sensitive to with with an IP experiment so this plot here what is showing you is a bit of a summary of what I just mentioned so here we've got each of the different geologic units mentioned so these are just the short names here so if we look for example the CLS is the Cambrian limestone and then units one two three and four are the mineralized units so the black curve here is resistivity so we can see that the cambrian limestone here is very resistive so that's that's sitting up here some of these other units are a bit more conductive so where we've got the hanging wall sandstone for example this guy or the siltstone shale are both on the more conductive side of things so that it's conductive because again it's a low resistivity but then here if we look at our mineralized units they're sort of sort of in the middle so it's not necessarily that DC resistivity is going to be diagnostic for these units but they're gonna help us build up the geologic picture and then in this case we can see when we look at the charge ability which is the bars here and it's on this axis here these four mineralized units are much higher than the others and so this gives us some confidence that charge ability is indeed going to be something that is constructive and will hopefully help us be able to delineate these units okay so now we're gonna look at our data so there as we mentioned we can collect both DC and IP data from the same survey so this is one line that we're gonna focus in on and this was the line of data it's originally published in that paper so we're gonna be able to go ahead and compare the results that we generate with sippin I to results that were published about 20 years ago at this point so what we're showing here is we've got what's called a pseudo section of of data and so we'll walk through in the code how this is actually built up and how to sort of conceptualize where a datum is plotted given the source and receiver locations but basically it's a way to image and plot your data where we can say that data that are close to the top are going to be probably more sensitive to the near surface where's the data that are closer to the bottom are going to be a bit sensitive to deeper structures and so visually sort of looking at it vertically that's what you can understand what's going on there sorry I am going to close my window because there's a lawn mower outside [Music] sorry about that okay so the the top plot what we're looking at what we've plotted here is apparent resistivity and so this this is just one way of visualizing our data we've translated electric potentials into an apparent resistivity and again in the code we'll walk through what that means this is these are not true resistivity values so you shouldn't be interpreting this as geologic units but it's just a way to visualize and understand and also QC the data and then similarly here we can plot up the IP data I'm still looking at these we can see there's potentially there's something interesting going on here we've got some higher charge ability here so there's there's something interesting happening in our data all right so now we've got these data our goal is to invert them to now obtain physical property models giving us the resistivity using the DC data and the charge ability from the IP data so what are we up against when we talk about inverting data so we've got our data and hopefully some uncertainty estimates that gives us a sense of how much we trust those data we also have our governing equations so we know for DC resistivity we're starting with a static Maxwell's equations and then we might also know something about the geology of the area so we might know what rock units were expecting to see maybe we have some well logs that give us some information about the physical properties of the rocks that we're going to be encountering so these are all of our inputs and then what we want to get out at the end of the day is a model and so in order to be able to do that we need to both be able to Ford simulate our data so meaning that given given a model of the earth we should be able to predict what what data we should be measuring and so for this we're gonna need to be able to break up the earth so write down some sort of discretization and then go ahead and solve solve those equations and then in the inversion what we're gonna do is now we're gonna set up an optimization problem where we're gonna try and fit those data so we want the data that we're predicting to match the data that we observe and in order to have that problem be tractable we're gonna impose some some regularization and we'll walk through what that means I so cific is what we've developed as a mechanism of solving these types of problems and so right now we've got utilities and functionality for simulating and inverting quite a number of data types so gravity and magnetics DC resistivity IP electromagnetics as well as fluid flow and so all of these problems can fit into a framework like this so as I mentioned there's two pieces that we need is we need to be able to predict data so we're gonna break up the world into blocks like Lego blocks and simulate equations over them so in this case we simulate the DC resistivity equations which looks like this and I'll walk through these and then as I mentioned before we then want to run an optimization where we're now going to try and estimate a model so we set up an objective function where we have one part where we want to fit our data and then one part where we use a regularization where we can impose some geologic constraints and also perhaps ask for a smooth model so that we actually can can recover something that's geologically interpretable and then we might have things like bound constraints so if we know for example resistivity should always be greater than zero we can't have negative resistivities then we might want to impose that so for DC resistivity here's the equations we're gonna be solving um and so this derivation is just to get you connected with the equations we're working with so we start from the electrostatic Maxwell's equations where here the e is the electric field and because it's curl is zero that means we can represent it as the gradient of a potential and so here we're going to use the symbol Phi as the description for electric potentials J is our current density and so in this case it's this is basically telling us a bit about conservation of charge is that the divergence of J should be the negative of our source currents and so our source currents if they're point electrodes so if we just have two electrodes in the ground we can think of those basically as Delta functions we've got one place where we're injecting current and one put one negative electrode where the current returns to you and then finally we have a law that relates the current density which is a flux to the electric field through the resistivity or the conductivity and so by piecing all of these together we get our governing partial differential equation so oops so we take the divergence of one over Rho times the gradient of Phi is some sort of so that's what we solve in the DC resistivity equations for induced polarization we adopt a linear model for chargeability for IP and so what we're going to do is we're going to linearize about the resistivity or conductivity and so here we use this symbol ADA to refer to the charge ability and what its effect does is if we increase ADA we increase electrical resistivity and so what we can do if we linearize if we linearize the resistivity is we can then actually write out a linear model for what an IP datum looks like and this is true here for a parent chargeability data which is what we're going to be working with for the century deposit and so here that those sensitivities this linear portion is just the sensitivities from the DC problem so what that means is if we want to invert IP data we first have to invert DC data in order to get those sensitivities and then we can go ahead and invert or simulate IP data so here where our governing equations for IP are just a linear set of equations and so the linear linear equations are similar to what you would encounter in potential fields problems as well so if you want to work with grab D data or magnetics data in those cases we're often working with a linear set of equations - so as I mentioned if we're going to work with IP data this is a two-step process so step one we invert the DC data to get out a resistivity model um that model and those sensitivities are then used to do step two which is the inversion of the IP data so we use those sensitivities then to perform this inversion and get a chargeability model so this is what we're going to be doing in in the first notebook and these results here that we're looking at are from the UBC inversion that was running 20 years ago so let's actually walk through what how are those results originally obtained so the image that you will see it on the right is published in the money paper so this was the original inversion an original image of the IP data on the line that we'll be working with so you'll see this number quite a bit as we go through the notebooks so that's that's the line that we're gonna work with there's 27 sources on this line and 115 or sorry 151 data in total these data were originally inverted with the dcpip inversion codes developed by UBC jiff in 1992 and the mesh that was used to run this 2d inversion had about a thousand cells and at the time it was run on a large Sun workstation it took about an hour per iteration so in order to obtain this image here that was about a day's worth of computation and so we're gonna see what we can do 20 years later and in Python well before we get there let's take a look and actually see what we can understand from the geology if we're if we're looking at these inversion results okay so the first piece we're going to look at is the DC inversion result so here what we have plotted what the color indicates is resistivity so Reds here are more conductive so it's a lower resistivity and blues are more resistive and so if we want to try and interpret interpret this image in the context of the geology that we're seeing how we relate that you have physical image to the geology through the physical properties so here we're looking at resistivity so let's come back to the resistivity measurements of the geologic units at Century so in this case we see we've got some resistors so we got some near-surface resistors here and out here and looking at our resistivity plot and physical property values it's quite likely especially at the near surface that these resistive units that we're seeing are these carbonaceous or sorry that Cambrian limestone so you can see here that that tracks out tracks out near the top and then as we have this fault here it looks like this unit does go deeper so this resistivity seems to be this resistivity image seems to be agreeing fairly well with what we expect to see for the cambrian limestone then we can also see a couple conductive features so we've got a bit of a conductor here and maybe something up here and if we look both at the the physical property values as well as our geologic section it would it would make sense if these correlate then with what we're seeing for the hanging wall sandstone so this guy here which is uh Hugh on this side and then if we look in the the center region it's in corresponding to this sedimentary um segment here and so these seem to be aligning fairly well so these lines are showing showing these faults there too so so that gives us some confidence that we can interpret this image in terms of the geology it's also encouraging to see that here where we've got this conductive unit there is a bit of a disruption and a discontinuity where we see the fault and similarly here this is this resistor that's starting to come in doesn't necessarily exist on either side so we're seeing some structural information from the resistivity image as well as potentially being able to understand where some of these geologic units are all right so now to the chargeability so we use that resistivity model now to run the IP inversion so we're working with the same line of data but now we're looking at chargeability so we have loop being low charge abilities that means that these rocks tend not to support a buildup of charges and then we've got higher charge ability is red and so those units tend to support the buildup of charges and so what's really encouraging to see here is where we're expecting to see these mineralized zones which are sketched out here and here as we are seeing a strong chargeable response and so these and these align well with what is observed from drilling and so this image where we're seeing basically being able to image the extent of the mineralized region this was an extremely compelling result because instead of having to drill and continually use a drilling program to delineate the extent of mineralization instead of that being the only tool that you have we now actually have the image from data completely collected at the surface that we can use to try and estimate the extent of extent of the mineralized zones so this was a real early success story of how images obtained from geophysics can really aid the understanding of geologic structures and a decision making in mineral exploration so in summary we started from dcpip data I collected over century we looked at the UBC inversion results using the DC code we were able to understand perhaps where some of these sedimentary units lie and so get some geologic understanding as well as observing discontinuities and structures that align well with the faults that we see this information was then used to obtain an IP to run the IP inversion and here we're able to image the extent of the mineralized zones from from the IP data so these are the data that we'll be working with in the first in the first notebook that we're going to open up um so our goal is gonna be to invert the dcpip data from center using simple and so the first question we want to ask is can we can we replicate these results from 20 years ago using a different code and perhaps our own different choices and parameters in the next notebook we'll go through and we'll explore some of the fundamentals of the physics of DC resistivity and this will hopefully give you some ideas about how you can set up your own forward modeling with sim pic and also explore aspects of the physics like the currents charges electric fields and potentials and so if you're looking at survey design type problems hopefully what you learn here will extend to two other methods that you might be interested in and then in the last notebook we're gonna go ahead and we'll look at some of the important choices in the inversion there's lots of knobs and parameters that you can tune and so we built up some some ways to go ahead and explore those and hopefully build some intuition of what they all do one thing I'll note is that in this case with the the layout of the notebooks there's a lot of code so we try to be very explicit and walk through step by step some of the things we're gonna probably gloss over a bit today but will be very active in slack and happy to answer questions my goal in doing this is really to first give some ideas as to really how this impact framework works and so instead of hiding things in utility functions and things like that we're gonna really work through all of the details so but if you have your own data file that you want to get loaded into something and we maybe don't have a utility for that walking through this will hopefully give you the the conceptual model of how you would do that for your own data here's a bit of an optimistic schedule potentially so I say that because I was only planning for four twenty minutes for this talk and we're already at 30 but that's alright that's how tutorials go so we'll set it a schedule and we will see if we stick to it but either way here's the plan so we'll start with notebook one and we'll invert the DCN IP data from century and throughout we'll be highlighting aspects of this impact framework so this should give you an overall tour tour of how to use impact to invert field data will take a bit of a break after that and then we'll have some time for a bit of Q&A and so again please feel free to drop questions in slack I already see that there's there's quite a bit of activity there I know some folks are potentially running into some challenges getting up and running on the the code and so if you do run into any issues today and want to just follow along feel free to also use the binder link so that that should work and run online and that'll allow you to follow along in real time and then afterwards we can always give you a hand with getting simpe and the notebooks up and running on your computer yeah so if you have questions both technical or conceptual feel free to please drop those in slack there's a number of folks who are helping out helping answer questions and they will flag questions and send those to me and we can answer them after the break then notebook - we're gonna be diving into forward simulation exploring some aspects of how to set up before it simulation and what sorts of choices matter in order to get accurate results and then finally a notebook three we will focus in on inversion and then ideally we will end with a little bit of Q&A I also wanted to include just a few links to some resources I know that we've got a broad audience here and so some folks might be you know DC resistivities in their wheelhouse and further folks these might be newer newer concepts and so here's a suite of resources on geophysics and inversions so we've got geophysics for practicing geoscientists this is a textbook we've used a number of undergraduate courses and is meant to be an introductory tour to a whole suite of geophysical methods EMG OSA is a bit more of an advanced resource on electromagnetic methods in geophysics and then with both of those we've got some lectures recorded and then specifically with respect to synthetic here's a few links for getting you up and running so we've got documentation discourse is a good place to go and ask questions as well as this in pink slack I wanted to give a big shout out to this impact team and specifically the folks who are jumping in today and helping answer questions here's everyone slack or yet the swung slack handles so if you're running into issues please feel free to tag any of us directly I also want to thank yog wali and Doug Oldenburg for both giving us getting a set up and running with these data and providing some context background on century and going through and giving us a review on these notebooks all right so before we transition I'll just keep this up for a moment so that folks can copy the link if you haven't already but we'll be starting here is where the github repository is for the content and then again as I mentioned slack is where you should come to ask questions all right so so what we're gonna do now is we're gonna go ahead and go over now to the github repository so I want to show you what that looks like and if folks have not already downloaded the code and want to this would be a good time so if you follow that link it will take you here to the simple organization to the transform 2028 tutorials so we've got the data here there's three notebooks and then the readme contains all of the information that you need so the live stream which is happening right now um so I guess you're already successful if you're watching this and then the link to the slides that I just used is right here and then the Conda environment you'll be using is this guy so if you haven't gotten set up and running go ahead and take a look here now at the usage section so in step one you need to get python hopefully you've already done that and then here we've got steps to get clone and create the Conda environment for this repository so please go ahead and follow these steps if you run into any issues ping on slack if you do friend of use and you want to follow along life I just use the binder link here and this will bring you to a web interface where you can run the notebooks okay so once you've gotten up and running we're gonna start by taking a look at the the first notebook so the inversion of DC and IP data at the century' deposit perfect so this is the line of data that we're working with so we saw this this mathil already um so this 46 800 East line is this is the one we're working with again here's our geologic section and the result there were aiming for so Ken can we replicate this results in and maybe update the color bar I can here so long folks cringing silently on the other side by looking at a jet color bar I've seen many great great posts on that alright so the first step we're gonna do is we're gonna import our dependencies I'm just gonna make this bigger to make sure folks can see it alright so there's I put in links here to some of the core tools that were using so if you're not familiar with numpy matplotlib or i pi widgets there's links to the documentation on each of these these are pretty standard Python Python libraries and then there's some things that were going to be importing they're specific to the simple ecosystem so this package called discretize what that's going to give us is all of the meshing capabilities as well as the differential operators so when we're gonna build up a partial differential equation discretize is what will give us the gradient divergence etc the pine hat solver library gives us wrappers on fast linear solvers so here we're gonna try and use Paradiso is our silver which is much faster than any of the the basic solvers like Lu and then soup egg is where we're gonna be spending most of our time and so sympathy original acronym stands for simulation and parameter estimation in geophysics so this is just core Python that we're importing um import discretize PI not solver and then within simpang so these now you're gonna see throughout and I'll explain step by step once we introduce them what each these modules do but we've got data we've got mappings things like the data misfit regularization optimization you'll see these come out come up throughout and most of this is related to the inversion machinery and then here now we're looking at the physics so we're gonna be working with so we've got both direct current resistivity which we import from the electromagnetics module its impact as well as induced polarization cool all right so we will import these guys so if folks are not familiar with Jupiter um what I'm gonna be doing when you see this transition from green to having a star in here and then down to the next cell that is I'm using shift enter on my keyboard and so if you're following along and want to be running the cells as we go that's that's what I am doing I'm so here we're just gonna set some bigger font sizes so we can see our plots alright so the first step we're gonna be doing is actually loading up the data so the survey the survey collected at century is what's called a dipole-dipole experiment so what this is we've got our dipole source so it's gonna be 100 meters apart we've got a positive electrode a and negative electrode knee and that's that's our source and there's a fixed seperation between those so these are gonna be 100 meters apart in in this survey and then we've got our dipole receivers so we call the positive 1 M and the negative 1 N and that is what we measure a potential difference over and so to probe different depths what we do is we're just going to move and separate separate both the source electrodes and the receiver electrodes by different distances and so basically the further further apart that they are the deeper that we should be imaging because we're looking at a larger volume of the earth so in this repository we've got a directory called century and so what I'm going to do is we're just gonna list that and see what's in there so there's a few things that we don't need to worry about dot files and stuff like that the ones that we're focused on for for this tutorial are the numbers and then where we've got a letter afterwards and each of those is now that contains the data for a specific line that was collected in the century deposit and the one we're gonna focus in on is this guy and so what I'm gonna do is now just print what's what's in that directory so there's a few different things in here I've got some log files some ODS files out files etc the ones that were most concerned about are the odds - those are our observations the rest of the files are actually related to the UBC inversion results so these are things that we're going to be able to compare to you so these are archived archive conversion results but the ones we're gonna be working with in this tutorial are this one so the P o T stands for potential so these are our DC resistivity data and then this guy here IP is our induced polarization data so these are the two data files that we are gonna be focused on okay so now I'm going to just show you let's just open up the file and let's let's see what's in the file so if we print out the contents of the file the first thing we see is we've got some information up top so this is just a comment it's a bit of a description it's a description of the data the next line here is the number of sources so we've got 27 sources in total and then both of these are just flags if it's a one it's a it's a dipole for the source and for the receiver so this is a dipole dipole survey alright now we're actually going to look and see what the data are so the way that this is structured is this first so the first line after we know how many sources we've gone is the location of the first source electrode so this is the location this is the a location so the location of the positive electrode source electrode and this is the B location so the location of the negative electrode and you can see that they're 100 meters apart this number here is then the number of receivers that are listening to that source or collecting data from that source and then so the next two lines are then the location of those receivers so this is the M location of receiver one and the end location of receiver one this is the measured potential and then this is the assigned uncertainty so these data the uncertainties have already been assigned for us in notebook three we'll explore the case where if you have to assign your own uncertainties but for this inversion we'll just work with the uncertainties that were used in the UBC inversion okay so this is what our file looks like so now we actually need to get it into a useful form in Python so what we're gonna do here is we're gonna write a simple function that is gonna read those data in and so within simpad we have lots of utilities to read in data files and things like that but here I wanted to do it from scratch so that you get a sense of how how we can actually work with these data so there's two things you can input so real input the file name and then verbose so in this case I'm going to print some stuff so that we can just kind of make sure that we've read things in properly and for now we're just gonna read it into a Python dictionary so we're not doing anything with simple yet we're just gonna read these data into a dictionary that's going to contain the locations of all of the electrons the observed data the standard deviations and the number of sources so here this is just a function in numpy that just reads reads the content in two strings so each line is then I know PI rate of the string with the like each new array being dilemma' dated by a new line so basically we're gonna end up with an array where we've got however many lines there are in this file that'll be that length of that array we know that the first line is the the first line after the comment line is the number of sources so we grab that if we're gonna just print some information we'll print that out and then we'll initialize some storage for the source electrodes and locations of the data so we know that there are 27 sources um so we'll just create that space directly for a and B and then the rest we're going to use lists because we don't know from starting at the very top of the file how many M locations for example that we have and then here I just have a little index where basically we're gonna just keep track of which line we've read so we'll loop over sources and then for each source we're gonna read in its receivers associated receivers and data and then here we'll just print out some information to make sure that we've we've read things incorrectly so that's all that this is doing is we initialize space for the receiver electrons that are associated with a given source and same thing for the data and then we just put all of that information into a dictionary so let's let's run that I said verbose is true so that we can kind of print some things out so here we ready and we've got the right number of sources for source 0 we got the a location B location and number of receivers so we can we can make sure and sort of double-check based on the data file that we are looking at if we've read things in yes I will make my text a bit bigger here perfect focus is that text size a bit better I was getting getting some feedback there okay so here let's now actually take a look at what some of the information is in that dictionary there we go okay so we've got now there are number of entries in our dictionary so our a location is an umpire as we said before our end locations are list saying the same thing as the end location standard deviations and observed data so let's just take a quick look at those so our a locations we got an array of all of our a locations so this is the the whole survey we've got it should be the length should be length 27 as we expected so then let's also take a look and get a sense of now what our receivers look like because these are a bit more involved we know we've got 27 sources in total but the number of receivers and the number of data associated with each source is a bit different so the way they were going to track them is now basically so we see that for example the end locations is a list so it's a list of numpy race and so the first array is the M locations associated with the first source and so here when we had two receivers for the first for the first source we should have two M locations associated with that first entry and this is the same organization we use for the observed data for example so same idea here so we get two datums associated with that first source so you can sort of see there's this hierarchy here not where we've got sources first and then under the sources we know how many receivers and how many data are a source are associated with each of those okay and so this is basically going to be mimicked and how we set up things in simple so so far we haven't actually touched simpe we've basically been in pure Python so now we get to this imp a part so the survey object in simpe is going to be what keeps track of our geometry of the sources and receivers in our survey so similar to the way we saw that that file was structured each source just gonna take your receiver list which knows about the different receivers that are listening to or associated with that source so just as a simple example here if we take that bat first the first two datums that we're going to get from our file we know we've got M locations so we have two and locations and two end locations and it's a dipole so we can create a dipole receiver so you'll notice here that I'm not creating two receivers I'm only creating one receiver for one receiver object and that's because it's all of the same type if you for example how to survey where maybe you had dipole and poles then in this case we would need to add the pole receivers to our receiver list but if you only have dipoles then we only need to create one receiver object associated with with those multiple locations another thing to note is that in simple we use column-major ordering so that's the same idea as Fortran ordering so here when we talk about the X locations those are gonna be in the first column now we're doing we're gonna be doing a 2d simulation and a 2d inversion so these need to be located in 2d space so we're not going to consider typography so these zeros here are the Z locations and those are just going to be on the surface of the earth if you did have to fog rafi for example then this is you would input those guys you didn't put the correct locations here okay so that gives us our first receiver and now we want to create a source so our source also has a location and now a source has a single location and so we're just gonna put in the a electrode so we know where it is a long line and then again we are not using we're not using topography so it'll be at the surface so we tell it where our a location is where a B location is and then we can go ahead and construct a source object and so it is a dipole source with an A location a B location and a receiver list and so very much similar to how we structured that dictionary that's how we're going to structure our survey in sympathy and then the survey object at the end of the day is then just a contains a source list so we've done this for the first source and then we would go ahead and repeat the process for all of the sources that we've got in our survey so now we're just going to implement that in code so we start we initialize an empty source list and now we're going to just loop over the number of sources that we had construct the receiver objects so see here we're gonna grab those m locations for each of the sources so this for the first entry will be length 2 for others it could be longer same thing for the end locations we'll construct our receiver object construct the source which has this location and the receiver list and then we'll have we'll return that source list so here if we then go ahead and look for example at the length of the source list it's 27 and if we look at a source object so let's just say we've got a source 1 is now source list at 0 we can look and we can see some of its properties so source 1 it has a current it's got a location and a location of the location a number of data associated with it so that should be 2 and then it also has its receiver list scoops it also has its receiver list which in this case is that that dipole receiver so then once we have this source list we can then go ahead and create the survey and that is then just a container for that entire survey geometry so if we look at some of the properties of the survey for example we can grab all of the a locations all the below Haitians we can do things like makes we shouldn't do that there sorry some of the API has changed so I'm trying to be careful what I draw your attention to so we can look at the number of data on the entire survey so that's 151 that's everything we need now to keep track of how were the data collected so that's just the geometry but the other thing that we need is we actually want to then keep track of the observed data and the standard deviations and so within simpe the way we do that is through this data object so a data object takes a survey which is that geometry we just created it takes observed data so in this case the ordering of the observed data is exactly that sort of same wrapping that we just did for the survey so the first datum should be associated with source 1 and receiver 1 the second one should be with source 1 and receiver 2 and then so on and so forth so it's the exact same wrapping so since this is in the same order as what how we built up that dictionary we can just stack the observed data and the standard deviations and everything is in the right order ok so now let's actually take a look at our data and so as I mentioned before we can plot data in pseudo sections so what I showed you in the slides was apparent resistivity but that is a derivation basically from the voltages so these are truly the voltages that are measured at those the potential differences that are measured by our receivers and so you can see it's kind of a little tricky to interpret we've got some dark values some large negative values you know at the surface and it gets smaller as we go deeper and so interpreting purely voltages or potential differences from from a pseudo section is a bit non-intuitive because what you're basically looking at is we're getting information both about geology so the physical properties of the subsurface as well as survey geometry because if you've got a source really close to a receiver that's going to give you a stronger potential difference whereas if you go further apart and those potential differences will be smaller so we sort of have all of this information conglomerated into one and so that's where we're gonna go to apparent resistivity but before we do that I actually want to walk you through how these data are plotted because when I actually first plotted these up this is a bit of an odd pseudo section if you're used to looking at DC resistivity data you often sort of see like a V because we have more sensitivity and more data near the sub near the surface and less spots as we go deeper just based on the geometry so I wasn't actually sure if this was correct or if I'd made a mistake somewhere and so I built up a little app to actually go ahead and visualize how we're building up the pseudo section so this is where I'm going to use the eye PI widgets library so I'm not necessarily gonna go through this code it's basically just a whole bunch of geometry of plotting where these sources and receiver locations are and I think it's easier to actually just see we visualize it so here the blue dots are the source location so the dots are actually the positive and negative source electrodes the X is the midpoint of those and then the orange ones are receiver locations and again the X's are the midpoints of those and then so when we go ahead and we plot plot a datum on a pseudo section what we're basically doing is we're drawing angles straight lines from those and and plotting that in in space we're pulling that at a given depth and so the deeper we are the further apart they are the shallower we are the the closer together the sources and receivers are and so this first row are basically the closest point where the sources the sources and receivers are closest together in the survey and then as we go further apart we go deeper so let me just move this along so you can get a sense so we've got source 0 that's what those data looks like if we move over we added some sources there so maybe when some of the processing was happening maybe we actually cut out some data from these points or from those regions so if we look here this number and spacing is how we're basically choosing what depth that's plotted at so this source and this receiver they are the closest we're gonna get to in in this given survey so the first source and the first receiver we plot by convention and is one and then as we go further and deeper the end spacing increases so it's not a true depth but it still is somewhat related to depth because as the source and receiver further apart we are seeing deeper and so we can't actually see when we go ahead and plot up all of these sources and receivers that it looks like our pseudo section is correct but there's perhaps some processing that removed some of the data that would would have given us that V shape that we're used to seeing so the other thing I mentioned is that you know looking just at potential differences isn't always intuitive and so looking at apparent resistivity can give us a better way to sort of interpret and understand what our data are and so what apparent resistivity basically does is we're trying to effectively remove some of the impacts of geometry and so the way that we do that is we say okay what won't resistivity would we expect if the world were a half space because we can compute that analytically and so this is the formula for that and so you can see that G takes into account the geometry and then Delta V are our data and so using this we then go ahead and we can compute the apparent resistivity which sort of accounts for some of that geometry and so note that what we're doing here is we're basically saying you know let's assume the world was a half space and see what resistivities we would get so you can't interpret these as true resistivities because the worlds not have space and there's variations in the physical properties and things like that so this is not a true image of resistivity of the subsurface there's still some elements of the physics in there but it does give us a way to QC the data see if there's any outliers and see if there's any interesting trends and so this is the plot that I showed you that I showed you earlier where we're looking at the apparent resistivity um a pseudo section so there's some regions so out here for example we might expect that should see more resistive materials when we when we run the inversion whereas in regions like this where it's a bit darker we're perhaps looking for something that might be a bit more conductive and then - I mean this is can be a good way to spot outliers if you see a single point for example that has a very anomalous value then then that might raise some suspicions that maybe there's something that went wrong with that measurement um another way to view data is to look at histograms and this can be really useful in general both for trying to assign uncertainties and estimate uncertainties but also in more complicated surveys so this was a nice 2d survey but if you go out and you start creating a larger 3d DC resistivity surveys it can be non-trivial to figure out how to actually plot your data and so in this case we can plot up a histogram and see a nice distribution of our apparent resistivity values so we're going to use this for is to get an estimate of what should we be using for a background resistivity when we are running the inversion is because that's something we need to we saw we need a starting resistivity where we start the inversion from and as well as a reference resistivity and I'll show you where that comes in so in this case I chose to use the the median you could equally use for example the log mean sorry let me add a cell so it doesn't jump on and we can see that those values aren't aren't drastically different but these are things that you can go ahead and explore another common choice would be to use the maximum likelihood so we can use a value in here and so all of these are valid things to look at and so some of that really does come down to your geologic understanding of the setting and yeah and so you can try those out but in this case we're seeing that you know a difference of about 1000 meters or so shouldn't shouldn't drastically influence our results okay so now we've got oh I see there's a question coming it might look timely so with the contour F implemented within this function how do you change the noise levels ie increase that's a good question so in this case we are so you could use the pseudo section to get a sense of noise levels and things like that in your in your data um often it's helpful to look at histograms or in some cases folks do you repeat surveys where we switch out what the sorry reverse the source electrodes with the receiver electrons and the measurements that you should get should be the same by reciprocity so those types of things can give you a sense of what what your noise levels are and we'll explore that a bit more in the in the second notebook so feel free to come back with another question if that doesn't answer your way what you were looking for okay so we've got our data we've got an estimate of what the background resistivity looks like and so the first thing that we're gonna wear the next thing that we're gonna do is now set up the forward simulation machinery so we know what our survey looks like we know roughly or we've got a guess of where we should be starting our model from so now we need to simulate some physics so these are the equations that I showed you in the slides and so this is where discretize comes in is we've got this PE that we want to solve and so in order to do that what we're gonna do is we're gonna discretize the world break it up into cubes and assign basically locations for where each of these parameters lives on that mesh so we've gone I've linked a couple another tutorial that Rowan and I wrote a while back that goes through this the derivation in a fair amount of detail as to how to actually derive the discrete version of the DC resistivity equations I'm not gonna go into detail here but I'm just gonna point out weird the different things live because this is important for understanding I'm just physically where where things put on a mesh so in general and basically everywhere throughout simple physical properties are discretized ad cell centers so you can think of basically the physical property filling the whole space in that cube or in for a 2d section in that square in this case we're gonna use a nodal discretization of the DC resistivity equations so we're gonna put Phi on notes you can also derive and drive the DC equations with Phi at cell centers and that's what we do in this other tutorial and so if you're curious you can play around with both of these they are both implemented in simple but for this specific discretization we put Phi on nodes so potentials are on nodes and then the electric field is then going to be on cell edges so between between nodes and so one important thing when you discretize and when we're actually solving when we're solving PDEs on a mesh is to know what our boundary conditions are and so in this case we use neumann boundary conditions which means that there's no flux of current through the boundary and so in order for that to happen what we need to do is make sure that when we design our mesh that the boundaries extends sufficiently far away that there's no current that should be going through that should be going through those boundaries and in notebook too we have a little app developed so that you can actually play around with some mesh design and see what happens when you either when you bring those cells or when you don't experiment when you don't extend the boundary far enough a couple other things to consider so I've got three three main items that we'll consider at this point and as I mentioned in the second notebook you'll get a chance to explore them in detail so the minimum cell size and so this basically controls the resolution and so what you need is when we're simulating these electrodes so we've got a source current going in one and current being picked up in another we need to have enough cells in between them to actually capture the physics of that current moving so if we only had one cell between them that might not be enough to actually capture the physical phenomenon and so we need to have potentially a few more and so as a rule of thumb we typically do about for at least four cells per minimum electrode spacing and again this is something you can play with in that next notebook so that gives us our minimum minimum cell size so that's decision one the next one is the extent of the core domain so basically how how big of a region in the mesh do we need to be using those small cells so the horizontal dimension is pretty straightforward we know where our where our electrodes are where the farthest apart electrodes are and so we usually use the extent of the survey and then perhaps add a few add a few cells on either end in the vertical direction what we want to do is make sure that the core region extends to depth which were sensitive to variations in the physical properties so as a rule of thumb this is typically a fraction of the maximum source receiver separation so maybe a factor of divide that by three should give you enough enough resolution at depth again that's that's something you can play around with and then finally the last item we want to consider it now is the extent of the full modeling domain so we've got our core domain which is really the region we care about imaging but then we're going to add some padding to make sure that we satisfy those boundary boundary conditions so in this example we'll use a tensor mesh sorry I did not finish this sentence in the in the notebook but basically what it is is you can define the you can find each dimension by just specifying the cell widths so a tensor in each dimension is enough to define the entire mesh if you're looking at larger domains we have the tree mesh implemented and so there's examples of that in the in the discretized documentation so if you're doing a large 3d inversion I'd recommend that you check that out okay so here's the decision one that we were carrying about is what's our minimum what's our minimum cell size so in this case we'll make it the same in X and Zed if you for example wanted higher resolution and Zed you can make these smaller but we'll keep it the same keep things simple so we look at the minimum separation between our source electrodes so we looked at the a locations and the beat locations and for this survey it's a hundred meters and we'll do four cells between there and so with that we're then gonna get a base mesh for the core region of the mesh with 25 from here by 25 meter cells so that's gonna be in the core region and then we'll find our core domain in the X direction by just looking at where where are our receivers and how far apart they're so that's that the extent the lateral extent of this survey now to find the the Zend extent of the court domain as i mentioned we want to look at the midpoints and find the the furthest midpoints so here i look at we find the midpoints at the a and B and then the eminent and then we find the separation and find the maximum of that and then here we choose to divide it by three so you can explore in the next notebook some some different different ratios there and see what happens so here we've got a maximum separation of 800 meters and so our core domain and Zed is gonna go from negative 266 to 0 note here that Z is positive up and that's basically consistent everywhere through simple so anything below the surface is always Zed Zed negative as I mention before sometimes we want to add a few extra cells on either side just to make sure that before we start doing in the padding or that I guess the padding cells don't to strongly influence the the physics and so here we're just gonna add a couple extra core cells onto that core domain and then we'll add some padding cells and this padding factor here so for electromagnetic and electrical problems we can actually add padding these are nice equations for that seismic for example you don't typically want to be extending expanding or changing your cell sizes because you're dealing with a wave equation I'm but here we're dealing with diffusive equations and so we can get away with expanding cells further away from our source and so here I'm going to add 10 padding cells in X and in Z so here as I mentioned for a tensor mesh we can define two tensors and defining basically the spacing in the X and in the Z directions and so this here is specific to simpe the way that we set this up so the the first part in X that we do so we define the cell sizing the number of padding cells so this is defining the padding parsh portion on the left-hand side of the mesh and then the negative is a trick to basically say we want them expanding in the in the opposite direction um so I'll show you actually just what this looks like so I'm going to show you a plot first and then we'll come back to what was just printed there so these are those padding cells so when I said the negative flips flips that is we start with we start with the widest cells at the very edge of the domain then we can see our core region here we got padding on this side and we've got padding padding too depth some of the things that we've printed so the total number of cells this is going to be a big controlling factor on how expensive your computation is so if this is hundreds of thousands of cells it's going to be very expensive on the order of thousands we should be able to run this pretty quickly some of the other things we're worried about as I mentioned is the extent of the mesh so making sure it goes sufficiently far beyond our survey the minimum cell sizes and then here we also show the maximum cell size so those are the biggest padding cells and so is this beautiful printing here this was a contribution by theater so if you were enjoying that I should give a shout out to him okay so now we've got where we're actually gonna run our computation so let's actually set up a setup and construct the the forward simulation so this is the simulation is kept track goes by a simulation object in something so it contains everything that we need in order to compute fields and fluxes on a mesh so we instantiate it with the mesh because that's our computational domain and that also gives us all of the differential operators the other thing that we're going to instantiate it with something called a mapping and so what this is in simple we've defined the Maps module which gives us a lot of flexibility for basically defining what our inversion model is so we know that the DC resistivity equations we need to be working with resistivity values on the mesh but that's not necessarily what we want to invert for so as I showed you earlier resistivity can vary on on several orders magnitude and so perhaps a more natural description is actually taking the logarithm of resistivity and inverting for the log resistivity and so if we want to do that then our inversion model are gonna have values of log resistivity and in order to get that to be something meaningful on our computational map in our computational mesh we need to take the exponential of it and so that's what an XPath does and so in the inversion it's important that we have derivatives of all of these operations because that's gonna be how we compute updates and so the mappings module keeps track of both how do we translate from R the set of parameters that we invert for to a meaningful physical property on the mesh as well as all of the derivatives and so once you start to get familiar with this mapping concept it can actually be pretty powerful so for example we use this to also make sure in electromagnetic inversions that we don't invert for air cells because we need to simulate the equations in the air but we know we know where the topography is and we don't want to be inverting in here so we can choose the subset of cells that we care about in the inversion you could also do things like defining parametric models and so if I wanted to define a layered earth model and instead of inverting for resistivities in boxes if I wanted to invert for a layer depths and resistivity of each of those layers we could construct mappings to do that but for now we'll just be working with we'll be working with the Rho map okay so the other things that we need are a solver and so in this case we're gonna use Paradiso which is a bit of a faster solver for symmetric systems and then the other thing that we're gonna do is set this storage a to true so for smaller problems where we don't have where we're not memory limited the inversion is gonna be a bit faster if we store the sensitivity matrix rather than computing it on the fly so for larger problems you don't necessarily want to store the sensitivity because it is the size it's a large dense matrix that is the number of model parameters by the number of data so if you have lots of data and you have a big mesh this can get very very big and so you may not actually be able to hold it in memory and if if not that's fine because we have everything implemented to do matrix vector but if we aren't memory limited in this case we're not it makes things faster cool okay so here as I mentioned we're going to use this X map so that we invert for log resistivity so we provide the mesh the mapping the solver the survey and storage a and in this case we're going to use that the nodal discretization there's also a cell centered discretization within simming Convention is always so a simulation and then the dimension so this case we're doing a 2d a 2d for it simulation and then we're where things are discretized so in this case on nodes so I just want to show you some of the things that we have access to once we've created a simulation object so we can do things like predict data so deep red is gonna predict data we can get the fields and so this is something kind of cool and we'll explore that in the second notebook where you can access the things that we solved for so you can get potentials everywhere on the mesh you can get electric fields and currents and charges so we can actually start to visualize the physics using this Fields object we also have things like getting the derivatives so computing the sensitivities and sensitivity vector products and such and then you can see that that row map that we set earlier as well as the solver so yeah and there's a whole bunch of other things that we use when we actually sold but the deep red which predicts data the fields and then the make synthetic data methods are probably the ones you'd want to interact with most as a user okay so now we've but we've set up everything and we should be able to go ahead now and run a forward simulation so the first thing that we're gonna do is we're just gonna run it forward simulation over half space and this can be a good way to go ahead and check your the check your mesh because if we compute pre compute apparent resistivities over half space though be the same as the true resistivities because that equation is valid for a half space so that's what we're gonna do here so we set up our model you'll notice it's log resistivity and it should be the same size as the number of cells in the mesh because we discretized resistivity on cell centers so the vector that we input for forward simulations and versions etc needs to be the same size as the number of of mesh cells in this case so we go ahead and run that and now we can plot up the apparent resistivities so here I've set this color bar limits to be what we used previously and so we can see it looks like a half space which is as good that's encouraging if we comment that out we can see there's there is a little bit of variability but our true resistivity the resistivity value that we input was 135 I believe let me check let's double check that again so is Rho not yeah so 135 136 so we can see is centered around there and we've got a variation of maybe 400 meters in each direction one thing that is is insightful to look at so we see yes there's some variability in the pseudo section but then let's actually just look at the error in apparent resistivity so I've computed the air so we've got the apparent resistivity from the data that we just computed and I'm just going to compare that to the expected half-space resistivity I multiplied by 100 so this is a person and so in this case we're within like two and a half percent for these which is it's sufficient so the thing that you need to look for here is that we want to make sure that these errors are less than the errors that we set in the data and so typically we're looking at errors on the order of at least five percent in the data or uncertainties in the data that are at least five percent so as long as we are sufficiently below that below that value then we should be in good shape because the errors on the on the data are larger and we're not trying to fit the data to within two percent for example okay so that's everything now we've got everything set up to basically go ahead on given a model we can now predict data so that's everything we need for the forward simulation so next we're going to go ahead and actually set up an inversion and so as I showed this I talked about earlier we're gonna set this up as an optimization problem consisting of a data misfit term which we're going to call Phi D and then a regularization term Phi M that is weighted with a trade-off parameter beta and so beta controls basically how important are the data misfit versus the regularization in the inversion now we're gonna try and fit data we want to reach what we call the target misfit and this is controlled by the number of data so if we've set appropriate uncertainties then we should then what we should be able to do is reach our target misfit which is n divided by 2 because we define Phi D as 1/2 of a weighted difference between our Ford simulated data and our observed data so here the way the uncertainties come in or the standard deviations come in or through this WD matrix so we're just gonna define this data misfit term so we've got our WD matrix and it's a diagonal matrix where we put 1 over the standard deviations along those diagonals and so what this does is it then says ok if I've got a really small standard deviation then this is gonna be a big value and so that means that if this difference the difference between the simulated data and the observed data if that's large it's gonna get amplified because we multiply it by something big and so then the inversion is gonna then try and minimize that we want to try and reduce that so that's that's the influence that that has if epsilon is really large then the opposite happens is that we basically reduce the influence of this difference and so if you know that some regions or some data are noisy then you can increase this under deviations and that'll reduce their influence in the inversion so looking at this equation what do we actually need in simpe to work with this so this is all through thou the data misfit module and we'll use an L to datum s fix we're doing an l2 norm and then it's instantiated with data so the DC data that's gonna be what gives us our diab as well as the uncertainties so the data object keeps track of both of those and as well as a simulation so a simulation then is what is going to be able to give us this term here so we go ahead and create that okay the next thing we're worried about is on regularization and so this is important because we're dealing with a problem that's ill posed an easy way to sort of see this is if we think about a linear problem and Matt wrote a really nice tutorial which I put a link here to is that this this matrix G has many more columns than rows and so what that means is we got many more model parameters than we do data and this is true basics I mean this is true always because even if we tried a parametric disk or disk realization description of M what we're doing is we're imposing a different set of assumptions and we're imposing sort of a different style of regularization as we're basically saying the earth´s can only look like this but in general we have the the general thing that we have is a continuous earth and so no matter how we break that up this is always going to be an ill-posed problem and so we have to impose some sort of assumptions to get back a reasonable result and so in this case the most standard and sort of approach that we're we're not imposing very strict assumptions these are fairly standard assumptions is using what we call it taken off regularization so it's composed of three terms for this 2d problem if we had a 3d problem we would add one more of these and so this first term is called the smallness term so we have this matrix WS which we can use to weight wheat the influence of these different model parameters and we take a difference between em our model and then some reference model and the reference model is important because here if we are thinking for example about resistivity and we're inverting for log resistivity if you say oh I'm not going to use a reference model I'll just set that to zero what that actually means is that we are then actually trying to get em close to zero right we've set em ref to zero but what that actually means is we're trying to say bring em close to one ohm meter which is actually pretty conductive and probably not geologically reasonable for a lot of the settings that we work in so it's important to think through and and set this reference model to something that you expect to see in the geology of the region the yeren these next two terms are what we call smoothness terms and so these are gonna be looking at differences between model values and neighboring cells so basically what we're trying to penalize with both of these is large gradients we want to say let's get out let's get out a smooth smooth model and we can interpret from there now the Alpha values as what we saw when we combined the objective function of the data misfit and regularization with beta in between these alpha values weight the relative importance of each of these terms in the regularization in a notebook 3 you'll get a chance to play around with those so here's the way that we set this up is we use the regularization module in simple and we're going to use the T knopf regularization which is of this form so we need a mesh because that's gonna be what we use to constructs these differencing operators the the gradient operators so that we know basically which cells are next to each other and then we can set the alphas I'm so you'll see here I said Alpha X and alpha y equals 1 now notice because we're doing a 2d we're doing a 2d simulation and inversion in simple the of the dimensions is always XYZ and if we are only in two dimensions then we just take the first two of those which is x and y so conceptually in and all the code you saw me write above I used C because we're thinking about it with this but when you actually set stuff up in 2d and simple just keep in mind that you need to be interacting with the Y component as the second dimension the other thing you notice here is I use a different value for alpha s and that's because the values of Alpha X and alpha Z have a different dimension as compared to alpha s so for taking a spatial derivative on these guys there's a length scale that comes in and so this is this is one mechanism for trying to rebalance them and so again in on notebook 3 you can you can play around with all those parameters and you're actually you're welcome to play around with them here too this is a good example to get a sense of how these parameters influence the results but notebook 3 has a synthetic example so we know the true solution where is here here we don't okay so we've got basically the components of our objective function so now we're going to set up an optimization so here we're gonna use an exact cows Newton and we set a maximum number of iterations so you can set this as high as you want um but for you know exploring things we want to make sure that we get we get to the end the end of an optimization run within a finite number of iterations and then this controls that the iterative solver under the hood and again that's something you can play with later on ok so now that we've got the data misfit regularization and an optimization defined we basically have everything that we need to come up with what we call like a statement of an inverse problem is we want to minimize our data misfit and regularization using this optimization routine and so that's what the inverse problem module in simpe keeps track of and it has all the capabilities to basically go in and evaluate um evaluate the value of the total composite objective function as well as its derivatives okay so we've got a statement of an average problem and now the last thing we need to do is basically say how do we want to solve this and so that is where the inversion module of simpad comes in and so it takes the inverse problem that we just set up so we've got data misfit regularization and optimization all in here and then we take what's called a directives list and so at this point what directives are is they can orchestrate updates during the inversion so for example I talked about the value of beta which controls the the relative importance of Phi M and Phi D we might want to try and estimate them and so in order to do that what we do one common approach is to try and estimate what are the largest eigen values of the basically the operators controlling each of these so if we make a change which which of these terms is going to have a larger influence and so that's what this beta estimated by eigen value does the ratio here that we inflate allows you to basically dial up or dial down the importance of the regularization um based based on that estimate so here what I've basically said if we set the ratio to 10 well money 0 is sorry what we set it to is 1 then then I'm saying we should treat these is basically equal if we set it to 10 that would basically say let's try and have the influence of the regularization be 10 times 10 times larger than the data misfit initially the Buddhists schedule and this is a pretty this is a fairly common approach in especially nonlinear inversions is what we're gonna do is we're actually going to update the value of beta over the course of the inversion so rather than starting with a fixed beta and saying let's try and fit our data that way what we're going to do is we're gonna say okay let's start with a high influence on the regularization and at each iteration as we get to a try and fit a smooth model because that's what the regularization is trying to do we're gonna try and we'll dial that back slowly as we start to introduce structure in the and so dialing back basically just means we're going to reduce the value of beta and allow the data misfit to have more influence in the inversion so the cooling factor it's gonna reduce beta by a factor of four every two iterations is what that means and then finally we're gonna use a target directive and so this is just gonna stop the inversion when we've reached the target misfit if we didn't use this it would run for twenty iterations because that's what we set is the maximum number of iterations in the optimization alright so we've got everything let's run it so you'll see here the first thing we noticed that popped up is we set the reference model to and not so as I showed you when we're talking about the regularization earlier we want that reference model to be something geologically reasonable and if you don't set it by defaults simple we'll set it to your your initial model but you can set them to be different things if you set the MRF property on the regularization there's a few things you can see in the log so we can see the beta value the Phi D the Phi n you can see it only took us three iterations to to converge to a solution to reach the target misfit so let's take a look at our predicted and recovered data so I'm plotting two pseudo sections here so the first one we're gonna look at what are the true observed data and then the second plot we look at what are the predicted data from from the inversion and so overall that looks like a good fit I mean we're seeing similar features in both pseudo sections and so it looks like we're doing a good job fitting the data we can also look at the miss Finn so this is looking at a normalized misfit so what we do is we look at the difference between the predicted data and the observed data and then normalize by the standard deviations so if we fit all of the data within their standard deviation these values would be close to magnitude one and so in this case if we're seeing most most values somewhere on the order of one then we're doing a good job fitting our data we also see there isn't anything that there's there's not too much structure coming into this data misfit pseudo section so we're not sort of it doesn't look like we're missing any major features in in the data cool then here I showed you those fault lines earlier in the slides and so this is just going to load those up so that we can actually take a look and and include them in our plots so let's take a look at our inversion result and so you'll see here I make sure to map the inversion result right because M opt we inverted for log resistivity but now we want to look at resistivity and so you can go ahead and we can actually see we're recovering some conductive structures like what we saw previously a bit of a resistor here resistive unit here and so we can go ahead and load up the results from UBC from the UBC inversion um 20 years ago which took a day and we can compare those and so overall you can see we're seeing similar structures the UBC one was run on a coarser mesh so we've got 25 meter cells these are 50 meter cells but overall from the eyeball norm I think this is a pretty encouraging result is to be able to see we've got these resistive features we're getting our conductor here this this moderate resistor is showing up them both and so that's pretty incredible is that two inversions with two different codes twenty years apart we can we can get similar results and so I see we're approaching I guess an hour and a half um so I think this is probably a good place to pause for about ten minutes give everyone a break go top up your tea and coffee and we can we can regroup in ten minutes and maybe take some questions then all right thanks everyone let's get started again after our break so just a quick recap is we've now gone through and we've inverted the DC resistivity from Center we compared that to the UBC results and to be very self congratulatory these are excellent results and they there's some nice agreement between the two of them and also so he came deserves a shout out for working through and doing the original inversion for both the DC and the IP data for century with simpang and so big thanks to Sagi for for working through both of these [Music] so next we'll move on to IP so we'll go through the IP inversion and then we'll we'll pause and take time if there's questions that have surfaced before moving on to the next two notebooks which will probably end up doing a short tour around given our timing but that's alright they'll be there for you to play with afterwards okay so to remind folks what we're looking at we're now gonna be looking at chargeability so basically the ability for a material to act like a capacitor so in this case what we saw for a century is that we're expecting that the mineralize units should have a chargeable response and so that's what we're gonna be looking for so as I mentioned before we you are using a linear as model for IP where we're using the sensitivities from our DC inversion and we're gonna be incurred inverting for a parent we're going to be inverting a parent charge ability data so we're gonna play a very similar game as to what we did before where we loaded up our observed data file so instead of this being a P ot has now just been an IP but then other than that is the same and so we're using that same read DC IP data function that we wrote earlier I'm just gonna turn the verbosity off because we sort of now know what's going on with them so this is just going to be that Python dictionary where we've got our a locations B locations M locations and locations observed data and um standard deviations so the the survey setup is very much identical to what we did in DC the major thing that you will see that has changed is that instead of using the DC module we now use the IP module and in this case we flag that our data are apparent chargeability data because there's a few different ways to represent and quantify IP data but in this case we're working with the parent chargeability if you're curious and are looking for a more background on IP I recommend both the GPG which I provided a link to in the slides and I believe should be linked in this notebook as well as the MGS I okay so we set up our IP receivers set up our IP sources so I can the only thing you're saying that's changed is that DC now goes to an IP and then we go ahead and can set up our survey as well as our data and then similar to plotting the apparent resistivities that we saw in DC we can also plot a parent charge abilities and so you saw this plot in these slides and so we can see some interesting features in this pseudo section folks who are familiar with interpreting and working with both DC and IP data would recognize this this pant leg structure here which is often associated with some sort of compact target it also looks like we're seeing something potentially chargeable out here so there's a few interesting there's a few interesting features in these data and just mentioning sorry to remind folks we're where we saw this briefly was in the in the slide deck and so so it's just showing both the field data for for these Suda's for this line okay so we've got now our survey in our data set up so we're gonna define the IP simulation and so very similar to what we did in DC we are now just working with the IP module we're still going to be using a nodal discretization the main things you'll see that have changed or that are new is now we are inputting the resistivity model that we obtained from the DC resistivity inversion and that's because we need this to construct those sensitivities so we need this to construct basically what our forward operator is the other thing is is instead of working with a row map we now work with an eight a map because the physical property that we are working with is chargeability which we use the symbol EDA for okay so that gives us everything we need now to simulate IP data so now we can generate and predict IP data given a model so again we set up a data misfit term and a regularization term like then using a tea knopf inversion in this case we're not gonna use an XML because searchability tends to vary on a linear scale sorry I should have included that plot from our slides but we don't expect chargeability to vary by orders of magnitude so in this case I'm taking an exponential is not an appropriate representation of that physical property so in this case we're gonna end up using what we call an identity map and so that just takes the identity so it takes the values as is and uses those but because we know that chargeability should not be negative we're gonna use a slightly different optimization than what we did for DC so in this case so for DC we used a Gauss Newton and in exact cows Newton with no bounds in this case we are going to impose bound constraints so upper bound we don't we're not concerned about but the lower bound does matter so we've got a zero lower bound we should not obtain negative charge abilities and so this will project any values that potentially try and give us a negative negative step in the inversion back to zero and then as before we set up the inverse problem with the misfit the regularization and the optimization so in this case we're gonna set a starting model and so instead of looking at a parent chargeability values like what we did in DC we looked at apparent resistivity values for I'm trying to estimate that starting model what we do here is it's pretty common in IP inversions to basically assume let's assume that most of our background material and most of the rocks that were we're working with are not chargeable and so it should often be it should be small or close to zero in this case we don't want to send it exactly to zero because that makes the optimization a bit challenging because that's exactly our lower bound there's mechanisms if you do actually truly want to set this value to zero there's ways that we can do that I mean if folks are interested in that feel free to ping on slack there's there's folks who can answer your questions and point you to things that you would need to consider like appropriate preconditioners and such and that all comes in with the directives so it's an interesting point but it is something that that you can dive into in more detail in with questions but for now we're just going to pick something that's that's sufficiently small and we're gonna use the same directives as what we did before um we're gonna use the estimating a beta value we're going to use a beta schedule here we're gonna cool every every iteration because we're working with a linear problem so we can we can cool every every iteration because we've solved exactly at a given beta at each iteration okay so we can go ahead and run the inversion and so in this case we're seeing similar values similar printouts as to what we saw before the inversion converged in in four iterations and took seconds per iteration rather than hours per iteration has per twenty years ago so things things have come a long ways and then we can go ahead and actually look at the misfit so let's see if have we done a good job fitting our data so as as we did for DC we're looking at a normalized misfit so we've got our predicted data minus our observed data and then we normalize that by the standard deviations so we want most of the values to be in on the order of one they don't all need to be one but on the order of one um so again here we're seeing relatively small values there doesn't seem to be any major stress that we saw in that initial Summa sections that are still coming through here so that gives us a bit of confidence that we're we're doing a good job fitting our data so we can take a look at the results and so here it's quite different than what we saw in the DC inversion is we're getting too strong chargeable chargeable units both fairly well aligned with where we expected to see mineralization and again we have the UBC inversion results so we can go ahead and compare those and so overall again we're seeing very nice agreement between the two two different algorithms implemented in two different languages run 20 years apart is I think a pretty exciting result and so and again here we're doing a good job delineating the extent of the mineralization okay so with done that takes us to replicating results from from 20 years ago and so I want to assign you some some homework or some some ideas that you can all play with after after this tutorial so here we only worked with one line of data a century but there are six lines in total and so using the functions that we've built up here the the sort of file formats for those same for those lines are the same as what we're working with here the survey geometries are gonna be a little different but can you even birth them all and so if you take on this challenge please just send along images in slack and we'd be really curious to see what you folks get and maybe down the line we can actually compile compile these results and do a bit of a comparison I think it'd be an interesting study so yeah please feel feel free to play with these and and show us where you come up with so with that that takes us to the end of notebook one um we've inverted eye data from from century and our greeing with these with these results here so I want to pause for a moment and just see if there are any questions that come through on slack on that folks want to surface for us to discuss and while we're while we're doing that we can go ahead and open up notebook - and so notebook - is gonna be the DC forward simulation notebook and it should look should look like this here um so we'll just pause for a moment and I'm going to take a look at slack all right okay I see lots of discussion happening which is excellent but I'm not necessarily gonna gonna sort through all of those questions it seems like most most questions are being addressed by folks um so we can go ahead here and we'll get started with the DC forward simulation notebook so in this notebook what we're gonna do is we're gonna build up a synthetic model that uses the same survey geometry from that line in the century survey and explore some things so we're gonna explore some things like mesh design and also show you how to set up and visualize aspects of the physics and so if you're using this as a teaching tool for example if you want to explain some of the fundamentals of DC resistivity um oh I just see one one question popping in and so I I will ask this Usagi's posted this from the from the slack Channel so the question is is it fine for the cooling factor to go negative what would that actually mean so in this case you do not want a cooling factor to go negative in fact it doesn't make sense it doesn't work for beta to be negative we don't want to work with negative values when we're doing an optimization when we're trying to minimize something because what that it would actually mean is we're effectively then trying to we're then effectively trying to maximize a value of something so if we put a negative beta we would be actually trying to basically do the reverse of regularization which isn't necessarily going to give this positive results and a follow up to that is with respect to what happens if you were to increase if you were to increase beta and so in that case you could but what that would be doing is you would then be increasing the influence of the regularization over the data misfit and so effectively what that statement would be is you're then trying to effectively increase the assumptions that your of goes at the cost of respecting the data and so that's not necessarily something that you would want to do because we are the data are our source of truth our assumptions are things that we think we know about the setting but we want the data to be the source of truth but there are some settings where you might want to play around with with those types of things and you certainly could set beta 2 or that beta cooling factor to be less than 1 and see what happens because that would simulate that effect of increasing there's another question here that's come in so how flexible is simple if I wanted to find my own cost function for example in the inversion this is a great question and I think one of the things that I'm most excited about with seeing what folks can do with simple is that it is very flexible is that you can define we've set up the objective functions to be very pluggable and composable so for example you can define a data misfit term that's actually composed of two data misfit terms and that lets us do a joint version so if we wanted to do a joint version between gravity and magnetic data for example you would then construct a data misfit for gravity data one for the magnetic data and then add those together with potentially some weighting in between so tebow a stick has been doing some really cool work in his PhD thesis looking at joint inversions and so you should check out his thesis when it comes out which should be should be soon as well as some of the papers that he's written on then also if you wanted to define a different norm for example and so in the regularization maybe you want to try and define compact norms and so maybe something like and try and mimicking an l1 room so we want to try and promote compact structures that's something that you can do as well and so Dom Forney so along with that work and his thesis as well and so there's some papers out about about doing that and so those are those are two to look at and get some ideas but if you for example wanted to include a data misfit term that we don't have so maybe you want to use a different norm on the data misfit maybe we want to know one norm on the data misfit although we don't have that implemented as long as you implement the ability for that object to evaluate itself given a model and compute derivatives of that we need the first order derivative and the an approximate second order derivative then you're then you're good to go and that can be plugged in directly perfect okay and there's a few other questions I see coming in asking sort of about extensibility or different receiver types different data types so in general we tried to design a very modular framework and so as long as for example if you wanted to define a receiver type for apparent resistivity in a DCM version rather than rather than potentials that's something that you could do what you would need to do is then basically the thing that we simulate in the physics is electric potentials so then on your receiver what you would need to do is simulate or compute the compute the apparent resistivity and then also include derivatives so that we can go the derivatives of that operation so that we can run the optimization so in general we've tried to design things to be very modular and very extensible so that you can you can try out your new ideas and different data types and things like that for reference a good place to go to ask these sorts of questions and if it's something you are unsure of is plugged into simple or if you're not sure where to go to implement things if you go to a simple discourse so soup egg discourse dot group we have a community discourse channel and so here you can come and ask questions and I see there's a lot of questions we need to get through as a team so and some folks asking some interesting questions so feel free to come and post post your questions or ideas here and we can try and support you and users can try and support each other as you try and implement new ideas perfect okay so let's get started on notebook two so here we're gonna set up a synthetic example and simulate some physics so what we're gonna do the same imports is what we did before so nothing particularly nothing new here and again loading up the DC data this is the same function as you saw before if you were doing this in a project I would recommend you actually create like a Python file so this is something that you can import but here if you wanted a refresher or want it to go in and poke around and change change some of these values you could certainly do that so we read in those DC data and now we're gonna construct our source list again so now we've gone from Python where we butter our data in a dictionary and now we're going into simple in this case one thing I'm gonna change is rather than working in the coordinate system of the survey for the synthetic example let's just work in a local coordinate system so the center is zero and that'll just simplify things like setting up the model okay so we're gonna do now is we're gonna explore some of those factors that I mentioned with respect to mesh design so some of these things like looking at the number of cells per electrode spacing number of core cells how deep we extend that core domain the padding factor and the number of padding cells so this function basically just wraps what we did previously in Python code by setting up the mesh and just lets by putting in a function we can easily just go in and change and swap out these parameters so we can see that when we build that build mesh it's gonna return a mesh and then it returns those core domains and I do that because it's nice to be able to just plot the core domains um and we can see here that this is basically the same mesh which is what we notebook and the next thing I do is I'm going to find a couple functions so I'm gonna find one where we just run a forward simulation over half space so we can provide the mesh that we just created the survey a half space resist a beauty value and then this nky and so what this is is the DC problem that we're solving it's actually it's like a two and a half d simulation and so to get the other dimension the third dimension which allows us to actually truly simulate like point point sources in 2d we use digital filters and so the nky is the number of digital filters and so this basically it's a parameter we can increase and decrease if you decrease it we're going to decrease accuracy but also speed-up computation because each one of these is basically an independent computation that we then are gonna integrate over and then if you increase this we should increase accuracy so it's a parameter you can play around with and get a sense of here and then finally this is just plotting code and so this is a mechanism for looking at the apparent resistivities given predicted data and the half-space resistivity okay so if we stitch that together then now we can actually run basically a forward simulation in three lines of code which is nice so we can run the half-space we set a half space resistivity compute some predicted data and then plot apparent resistivities and so here this these are the two plots that we were looking at before this we were looking at the pseudo section where we know if we input a resistivity value of 100 we should be getting back basically a half-space response at about 100 ohm meters and then these are the error this is the error so this plot is looking at the apparent resistivity minus the true half-space resistivity divided by the half-spaces resistivity and then times 100 to give us that value in a percentage so by looking at this point we can then sort of mess around with simulation parameters and basically see what that does see what changes our accuracy so I stitched all of that together into a single function that basically takes all of those parameters so that we can build a little widget app around that so here this is basically just the code to define the widgets and give them nice names one thing you'll notice we use the continuous update is false because if you in if you input a slide bar for example each time you slide over it's gonna run a new computation and if we're I mean these are relatively cheap computations but they're not they're not completely instantaneous and so we'll each computation is run only when you sort of stop interacting ok so then before we start playing with this I'm just gonna pose a couple of questions and so let's look at those three decision factors that we had before is if we change the number of cells per electrode spacing padding and then the core domain extent okay so here we've got a few different slide bars so these ones change the cell size so this is right now we can interpret this is this is four cells per hundred meters in the x-direction so that's exactly what we had done before is we've got four cells per electrode I'm sorry per electrode spacing and then this is the same thing in the same direction these slide bars influence basically how many core cells do we add outside of the outside of the survey area this parameter here controls the depth of the depth extent of our core domain and so here this is basically one third of the maximum separation between our a our source and receiver electrodes this controls the padding factor so basically how quickly when we looked at the mesh sorry I'm gonna scroll up to the point basically how quickly those cells expand so in this case we used a factor of 1.3 but that's something that you can change this changes the number of padding cells this changes the half-space resistivity on a log scale and then this changes the number of digital filters so here let's go ahead and let's see what happens if we decrease the number of X number of cells per electrode spacing and so in this case we should expect that if we decrease this value by too much we should decrease accuracy because we're not able to simulate the physics as much and you can see that starting to break down so if we go only one cell per per electrode spacing um we're off by probably on the order of 3 4 percent or something um same thing now if we reduce the number in the z-direction we can see this was fairly dramatic if we go to course in the said direction we're now getting errors on the order of 8% and this is unacceptable for an inversion where we want to be fitting data to within 5% okay let's also take a look at the number of padding cells so oh sorry one other thing I wanted to point out let's go back here if I decrease this guy we see that for the most part it's really these short separations that are affected so it's the near surface where we're doing a very poor job because that's this is where when we have services and receivers close together that's where we're most sensitive to is the near surface and so if we're really doing a poor job sort of capturing the physics on these shorter length scales we'll see it in the short electrode separations whereas now if we start to mess with the padding which is going to control if we are satisfying our boundary conditions or not so let's go ahead and decrease this guy let's see what happens when we change the padding in the x-direction so you can see now the the short spacings aren't as affected but these further offsets now are and so that is that's because we're starting to violate boundary conditions is that we have said that no currents can go through the boundaries but the boundary doesn't extend far enough and so we're getting currents out there and so so the solution is breaking down and you can see even with seven padding cells so initially when I ran this I tried with six or so you can see that this this isn't this isn't quite gonna cut it if we want to get within five percent five percent accuracy when we're trying to fit our data so I encourage you to take a look at this and explore the rest of these parameters and and see see what kinds of meshes you can design and see if you can build up sort of your understanding of where where accuracy plays a role where these person area where these parameters play a role in accuracy the last one I'll point out is this nky so if we dial this down the computation is gonna get cheaper but it will likely get less accurate and so we can see where that starts to break down um and it's kind of happening all over all over the map and that's because this is really changing how we're solving the physics and so if we don't use enough filters we're not actually capturing the physics well enough whereas if we are using here we use we chose to use eleven and that's currently the the default we have fairly smooth errors with spacing if we increase this event we should be able to actually reduce these errors and increase accuracy so you can see that ax is changing and it's a touch deceiving when they're all updating but as we continue to increase this we should see a bit of accuracy increase and that seems like there's there some accuracy increase in in the middle regions here but perhaps not as much of these more extreme offsets so but between sort of eighteen and eleven we're noticing there's not a huge difference so overall we're doing a pretty good job with eleven and that will then speed up the computation because we don't need to do an evaluation for each of these filters okay so the next thing I want to show is actually going ahead and if we want to build up an example where we are exploring some forward simulations and want to explore the physics and so this has been a really powerful teaching tool and explaining basically what are the fundamentals of DC resistivity what sorts of things what sorts of geologic setting this does DC resistivity is it not conducive to you and what sort of targets are we very sensitive to so this we're going to just define a simple model oh we got a background we have a resistive block and a conductive block and then we just define the location of those so now what I'm going to do is I want to show you basically how you now put this geometry onto a mesh because now that's that's where we need it to be able to evaluate the the forwards emulation so I started creating a vector that's the same length as the number of cells because we need a single resistivity value for every cell and we'll just start by treating it all as the background and then we're gonna find the indices that correspond to each of the blocks so indiscreet eyes the way that we can find these is through the grid properties on a mesh so here you'll see some of the things that you have accent to when we have a mesh we have things like the area we've got operators like averaging we have operators like the cell gradient for example but the ones we're looking at here now are the grids so grid see see there's the grid for the cell center locations because also see we've got a X so that's the edge the X engines the faces tax faces and n is for nodes and so if you need to be trying to figure out where to put your physical properties or if you wanted to for example look at the electric field which we know is evaluated on X edges and we're trying to really compare do a code comparison or compare with an analytic we want to compare the solution exactly where it's evaluated so we might want to try and play around with the DX script to find the locations that we care about so again as I mentioned we use like the Fortran style of ordering so each column corresponds to a given dimension so in this case if we look at mesh grid TC dot shape see is 3900 which is the same as our number of cells by 2 and 2 is our number of dimensions so the first dimension here is X and then the second dimension here is what we're treating is it but again if we actually wanted to access that grid it would be called Y in discretize so here we find the X locations and then find the locations that are greater than the minimum and smaller than the maximum and then we do the same thing with zid and then we can define row at those indices to correspond to the resistive block and then we just play the same game with the conductive block so we can take a look at our model now very simple model so we've got our resistive block over here and our conductive block over here in this case we're gonna use an identity map we're not running an inversion so we don't really need to worry about translating between log resistivity and resistivity let's just work in resistivity space and then same very similar setup as before so we're gonna run two forward simulations we'll do the half space and then we will do the block model and this will allow us to basically compare how sensitive are we to the block so we can take differences within without the block so we can plot both pseudo sections notice that these are on very different scales apologies I should have put them on the same scale but this is basically identical to the background and then here we see sort of two of these like pant leg style um features in our pseudo section near the horizontal location where we put our blocks so we definitely do see an impact in the data so the next thing I want to do is we're gonna go ahead and we're gonna plot some of the physics and so you'll notice here sorry I'll scroll up when we actually computed the solution rather than just doing the deep red or the make synthetic data like what we did in the previous notebook here we actually access the fields object and so in simple a field object actually contains our solution everywhere on the mesh when you're running an inversion this is an intermediate step is we only need our solution at the receiver location is to base it to compare between our observed data and our predictive data but if we want to explore the physics and try and get a sense of where the currents are going where charges are being built up then we want to actually see that solution and have access to it everywhere on the mesh and so that's what the fields object does coming back to where things live this is just a quick refresher so physical properties live at cell centers potentials right knows electric fields and current density for this specific discretization live on edges and so for plotting purposes we want to average everything to cell centers so there's two averaging operators will use this have NTCC averages a scalar value from nodes to cell centers so we'll use that to a ver egde the potentials to cell centers and then e and j are vector quantities that live on edges and so we want to use something that preserves the vector quantity so here we use the average edge to cell center V and so the V preserves the the vector so we'll get three values and xly well sorry in a 3d mesh you would get three values and XY and Z in a 2d simulation we'll just get the two components so that's all this code is doing and then we we also include some plotting code so you'll see some things in here like I said different color bars it's nice to sort of see charges on red red blue but this isn't really anything too too involved it's just it's just plotting code so we'll go ahead and let's just take a look unless for example let's look at the secondary potential for source 15 and so you can see some interesting features actually let's maybe start soulless so if we look at the primary this is the background so just the half space model you can see most of the action is happening near the near the sources if we look at the total which is gonna be with the blocks don't see much different um but if we look now at the secondary it's the difference between those two and so rather than having to type this every time what we've done is I've just assembled this into a little widget but this is gonna be the only plot that updates and so we can look at the model on a linear or a log scale what I've been showing you before is on a log scale so that looks more familiar with the previous plot in the notebook and so let's go ahead and actually kind of just step through what's happening in a unit experiment so I guess this isn't showing up very well but the source location if i zoom way and sorry I'll fix this later on the source electrodes are here in here and you'll be able to see it when I turn on the current so here's our model let's go ahead and actually look at what happens when we turn on this source so here we're gonna look at the source electric our source current density this is on a log scale and so in this case the white space here this is a mistake I made in the plotting code above all also update that I should have extended it further but it's it's just an artifact of plotting so you can see here we can tell which one is the positive electrode and which one's a negative electrode so you can see the arrows are going this way so we've got our positive electrode and our negative electrode return electrode um you can sort of see it it looks like currents are being a diverted a bit around here and maybe perhaps preferentially flowing here you can see that Bend so that gives us an indication we got something conductive here and so what happens when we flow current and then it crosses boundaries between units of different conductivities we build up charges and so in this case I'm going to flip over to the charge density but before I do that I want you to build up a picture of what do you think should be happening um if we've got so for example we've got currents flowing around this guy currents going this way current preferentially flowing into this this target let's let's take a look um so I should do this on a linear scale stole my punchline and so we want to look at charge density on a linear scale and so basically right now we see that the the main charges in the system are due to the where we actually have our current electrodes it's a big plus right next to the positive electrode and a big minus near the near the negative electrode but if we actually want to look at the influence of the blocks we can look at the secondary charges and so we see some interesting buildup of charges on this was the resistive block and this is the conductive block and so one of the things that would be interesting and I challenge folks on slack to actually play with this model explore some of these physical physical responses and let's see if we can actually build up an understanding of what's going on so for example let's start here why do we see a positive charge buildup on this interface if we come back to the currents and now we're looking at let's look at total currents on a log scale it's a little easier to see what's going on okay so that the location of the block is right around here and so we see we have currents flowing from what is a conductive walk into a more resistive background and so basically we've got charges that are going through through a conductor they hit something resistive so you can almost think about them sort of like getting stuck there and so that is what gives us this positive charge build up on this interface here is when we cross from conductive to resistive R we get a positive charge both of them so one of the things that I think is constructive to your explore and really trying to understand the DC resistivity surveys and also some of its limitations is what happens if we add a layer about the blocks so let's say define a layer between said is minus 50 and minus 100 and start with something resistive where do you see the currents in charge building up do you think we still see the blocks in the data and then try with a conductive layer what changes and so if you want to post some of those images on slack then we'd be really curious to see see what you come up with and also fine come up with some creative models if you've got an interesting geology or an interesting setting that you're working in and you've got a conceptual model of that um see if you can plug it in and see what you see so feel free to send us images and screenshots on slack with your short explanations of what you think is going on and we're happy to chat through that with you and so here we'll pause for a moment before switching over to notebook number three and if there's any questions that folks have feel free to send those through and the simpad team on slack will will flag them for me and you know in a brief interlude while we wait for any questions that come through I want to point folks to the simple website so we've got the website is simple X Y Z this is our main landing page there's links with there's just a bit of a description of What's in page is walk through some of the components there's links for installing joining the discourse there's a video and so Rowan and Saki gave a talk at Syfy a few years ago that provides a nice overview and then lots of links for for getting connected the other links that we have here are for the documentation so the link is Doc's MPEG X Y Z and so here you can find some instructions for getting up and running there's a whole suite of tutorials so Devon Cowan deserves a lot of credit for bringing these tutorials and getting them into shape so that you can walk through and for example if you want to see how to forward simulate or invert magnetic data for example there's tutorials on tutorials on how to do that so referring back to that previous question about normos for example so here's one where we show a sparse norman version with with magnetic data um so you're welcome to check those out okay I don't see any questions coming through at this point so let's get going on notebook three okay and so sorry I scrolled right to the bottom of this one when I was checking the things we're working so we'll scroll up to the top I'm so in this one we're gonna use a synthetic example so the same or a similar block model as to what we used before to explore some aspects of inversion so it's for setting uncertainties and show you how to do that adjust some of the regularization parameters and look at trade off parameter look at setting a reference model and a lot of these choices and seeing what impact that they make in the inversion um okay so let's start with again the imports that we've done before and we'll go ahead and we'll load up the DC data so this is again the function that we saw before um again if you were actually doing this as a library you'd want this to be something that you store in a Python file um but for completeness and if you need a refresher it's all there and then we set up the survey again so we've done all this now a couple of times we build our mesh and so this is just the the same function that we used in notebook two and so we can go ahead and plot that and take a look so this is a mesh we've been using throughout and then we're gonna use the same the same model geometry that's what we did before so we've got our block model I changed sorry change the color scheme on you between notebook two and three just wanted this to be a bit more in line with what we were looking when we were looking at the century data okay so we start by defining our true model and in this case we're going back to Lord working in log resistivity because now we're running inversions again and we want to work with the log resistivity in the inversion so that it's strictly positive and it can vary logarithmically and that's the space in which we search for a solution so now we're gonna generate some synthetic data and so again we use this X map this nodal simulation and we're going to store the sensitivity so that the inversion is a bit faster now to generate the data we're gonna use this make synthetic data method and so what mate what you can do with make synthetic data is it's actually gonna return us a data object so that's the same thing that we were working with in the first inversion to capture everything about our field data so it's gonna have the survey in there the observed data as well as uncertainties and we'll assign those uncertainties afterwards in this case we're gonna add some noise to our data and so we're gonna add 5% noise in this case that's the default if you wanted to change that we can look at at making synthetic data so we can look at the documentation there and you can see that you would change the relative error term here in mean that makes synthetic data function so if we wanted to add 10% noise you would just add 0.1 okay so we generate our D yeah and then we can plot a pseudo section of those data so it's similar to what we saw in the previous notebook but there is a bit more noise and jitter to it and that's because we haven't done so as a quick refresher this is what we're doing to solve the inverse problem is we're setting up an optimization problem how we define a data misfit a regularization and a trade-off parameter between the two we're gonna run the inversion until we fit until we hit this target misfit and then potentially include some bound constraints in this case we're not including those but for the general case we could and then again here's what those terms look like is we've got our data misfit turn where the uncertainties are captured in this WD matrix and then our regularization term where these WX & wz are basically approximating derivatives is we want to enforce now there that there are not large gradients between adjacent model cells and the WS in a typical case is it could be an identity matrix but also if you are doing for example a potential fields inversion you could use depth weighting there or sensitivity weighting so there's a bit of flexibility with what you do with this and then the parameters that will also focus in on are these alpha s alpha X and alpha Z which influence how important each of these smallness and smoothness terms are in the inversion so the first thing we're going to look at is focusing in on the data misfit term and assigning uncertainties in the first notebook they were given to us but that's not always the case and so in this case what you want to do is we're gonna start we're gonna look at the histogram of our synthetic data typically when we assign noise or assigns three standard deviations we think about it in two parts so we think about it having a relative error so that's like a percent percent error and so in this case that's going to multiply the absolute value of our observed data and then we talk about also a noise floor and so the noise floor effectively quantifies what is close enough to zero and so it's important to include this when we have data that has both positive and negative values which we do in this case in DC we do because we don't necessarily want to fit you know if we have a small datum if it's close to zero we're not necessarily concerned about fitting that to ten orders of magnitude so the noise floor is what basically says this value is close enough to zero and so a good way to try and assign that is by looking at a histogram of your data because that gives you a sense of what is small so in this case if we look at the the log the log of the absolute value of our observed data we see a distribution like this so we've got basically three and a half or sorry two and a half orders of magnitude between our smallest values and looking at this so we see our smallest values around a log of 10 to the negative 3.5 so something like 10 to the negative 4 is small for these data and we'll ask us to fit all of these data if for example you didn't necessarily trust your smallest potential values which is sometimes the case sometimes we don't necessarily trust that really the really small values because they're very long offsets and we're not necessarily sure the instruments are always sensitive to them then you might want to define a slightly larger floor and then we wouldn't we wouldn't be fitting those data and then there's also a relative error so sort of how much background noise do you assume 5% is what we'll do in this case but that's again something that you can you can play around with and so just to show you basically where our noise floor lands is here so that is fairly small compared to our data so now that we've made these decisions so we've decided what our relative error in noise floor we want to assign those to the data so we can assign those by setting the relative error and the noise floor properties of our data object and then what I want to show you is that actually so within simpe what we use to construct the WWD matrix is the standard deviation vector and that's constructed exactly as I showed you before where we've got the relative air times the absolute value of the observed data plus a noise floor and so this little statement here is just to show you that that is the case is that these are all close is if we evaluate the uncertainties like this that's exactly what is going on when we ask for the standard deviation from the synthetic data object this next button um it's not necessarily super intuitive but I just wanted to show you to give you a sense of how the percent error and noise floor influence different parts of the data so what I've done is sorted all the data by magnitude so the the x-axis here isn't super meaningful other than larger numbers are bigger amplitude data than smaller numbers and then this blue shadow is basically what the what the uncertainty envelope is around those data and so we here when we're looking at large data it's the percent error that are the percent standard deviation that controls how wide that is whereas when we look at smaller data it's the noise floor so the 10 to the negative 4 that controls basically how large that's how large that uncertainty is so just giving a sense of based on the magnitude of the data which of those terms is is influencing things all right so next I'm going to create a couple of functions that set up the inverse problem for us in a way where it's gonna be easy to go ahead and change all of these parameters so this is very similar to what we did in the first notebook is that this function here is gonna be responsible for constructing the inverse problem and so that's the object that contains the data misfit the regularization and the optimization the parameters that can go into that that we might want to change so on the data misfit we can mess with the noise floor and the relative error so changing our standard deviations in the regularization we can change those alpha parameters as well as the reference model and then in the optimization there's two parameters we can change the maximum number of iterations and then the maximum number of solves and this can basically act like a bit of a regularization term as well as the CG solves are somewhat smooth so if we do less CG solves the results that we get are a bit smoother so that's something you can play around with it it won't be the focus of our next examples so here we set up the data misfit so if we've changed those relative area noise then we update them here and then we assemble that data misfit same thing with the regularization this is what you saw before the only thing to know here is that our reference model needs to be the same size as our is our number of cells in the mesh because that's the number of parameters that we are inverting for so in this case I want to perhaps just change a scalar value for the reference model but when I input that and instantiate the regularization it needs to be a vector that is the same size as the number of model parameters and then we do an an optimization set that up and assemble in the in the inverse problem one thing I wanted to show you and this is a bit more of an advanced step and so I'm not going to spend time dwelling on it but I wanted to give you a sense of how the directives actually operate and so here we're going to do is I'm actually creating a new directive class and so when folks were asking about you know sort of extending Simpa and exploring your own ideas this is one example of them and so here I inherit the base directive and then there's two methods that I'm going to cream so what I want to do in this one is I want to just save information that I might want to plot after the inversion runs and so what we're gonna do is there's two methods we define there's the initialize method and this is called when we first start running an inversion so when we did that him dot run of a model this is called before we do any sort of computation about what predicted they are or anything like that and so for example if in the like estimate beta directive this is where we estimate an initial beta and assign that in medium burst problem the end inner is run at the end of every iteration and so for example the target misfit is evaluated them so once we run to the end of an iteration the target misfit will just check and say like hey did we reach our target misfit if yes let's stop the inversion if not keep going the beta cooling is also evaluated at the end of an iteration and so um it's gonna check should I be cooling beta at this iteration if so it will reduce the value of beta and then you keep going in this one when I initialize I just create an empty dictionary where we're gonna just store some interesting inversion results so we'll store a thing like the iteration maybe the Phi M Phi D values the predicted data at the model etc and at the end of each iteration I'm just gonna append those values so that's all this code is doing is this just saying look at the current value and let's store it in this dictionary so there are utilities for doing these sorts of things in symphony and also for saving it to output files and such and so in general before you start implementing your own sort of custom custom piece that's worth taking a look at what's there but just wanted to show in case you do want to I'm so in addition to the saving results directive we're gonna use the the three others that we used before so we'll estimate an initial beta we use a beta schedule to cool the value of beta and then we will use a target missed fin and so one thing to notice with the target misfit and we didn't talk about this in the in the previous notebook but I'll just highlight it here um so one of the things you can do so we know a Phi D star should be the number data divided by two but for example you can include a chi factor that scales them and so if you multiplied this by a factor of two then we would stop when Phi D equals n or if you reduce this then we would try and fit the data when you try and reduce the value of Phi D more but by default that's that's one but I introduced that here as something that you can play around with if you want to explore them so this this function is going to take an inverse problem and set up an inversion so there's parameters for the estimating the initial beta there's the option to cool beta or not so if you want to play around with a fixed beta you can do that we can change the beta cooling factor and the rate so right now what this says is we're gonna cool beta will reduce beta by a factor of two so if beta was ten at the first iteration at the end of the second iteration um what did I say it was ten so fader was ten at the first iteration at the end of that iteration we would reduce it to five and then finally we said are we gonna use the target misfit as a stopping criteria or not and in some of the first examples I'm gonna show you we're not gonna use that because I want to show you what happens when we over fit data so this function is just doing exactly what we've done before we set up those directives we assemble a directive list if we're gonna use the target we add that to the list if we're gonna cool beta we add it to the list and now we then just have an inversion and I also asked it to return the target and the save because we're gonna pull parameters and extract parameters from those afterwards okay so let's go ahead and run run some inversions just given track time okay so we'll start with a row naught that is equal to the background so that's our starting model and then in this case we set up those we set up those functions for now I'm just using the default values like what we set earlier for both and we're gonna explore some of the things that happen here so I set my initial model remember we're inverting for log resistivity so I'm going to take the log and then multiply by the number of cells cuz we're in braiding in every cell and then we'll run that inversion and then here this is new is that that inversion log which is that that save inversion directive so here when we created the inversion we got back an inversion object the target misfit object and the inversion log which contains that dictionary of our results okay so we'll go ahead and run that and you'll notice in this case I said let's not use the target misfit so we're gonna run this in version four twenty iterations so it's just going to chug along for twenty iterations irrespective of if we reach the target misfit or not um so some things you can see is we can see that we started with a large beta here and it's decreasing every iteration because that's what we told it to do by a factor of two 5d is decreasing with iteration that's a good thing that means we're doing a better job at each iteration of fitting our data and so you can see you know we started at ten ten to the three and we're now reduced down to something that is on the order of one and if I am does the opposite fine start slow and increases because we're introducing structure and by introducing structure that means we are introducing gradients and we're introducing values that are different than the reference model so now what I'm going to do now that we finished our inversion run for twenty iterations I'm just defining a few functions that are going to be useful for plotting and analyzing the results so the normalized misfit like what we looked at before gives us a sense of how well we're fitting our data and you can see these values are pretty small because we ran for twenty iterations so this is smaller than what we've seen perhaps in other inversions looks like we've gone past the target misfit then the next thing I'm gonna do is now define some functions where we're actually looking at those inversion results and we're gonna plot things like Phi D Phi M etc and so in this case we can plot the inversion results um either as a function of iteration or as a function of beta and look at those we can plot the misfit and regularization and we can plot the regularization components so the smallness smoothness and such as well as the TK off curve so rather than explaining all Cote I'm just gonna show you the images and then we can revisit them so we have these functions and so the first thing I do I created three plots so in the first one we will plot beta as a function of iteration the second one will plot the misfit and regularization and then the third one will plot the regularization components okay so we can see we started off on with a large beta and we reduced beta by a factor of two every iteration so that's what we can see in this plan this pot here is showing us now on Phi D which is an orange so our data misfit Phi M which is in green you see those are on different axes so Phi D is on this axis here if I am is on this axis here because they're on different scales and they represent different things but beta is again what relates them and then here the black line is the target misfit I also included a black dash dot line if you want to play around with this Chi factor and if you do that you'll see that these lines separate and go to different places but for now we're treating them as the same so you can see after what is it one two three four five six sessions it looks like we fit our data um but here we continue on and we can look at those model results and then here this is looking at those the smallness and the smoothness terms and seeing seeing their influence and I've included the Alpha values in that evaluation um so you can see here it looks like the smallness term increased the most and so we're getting structure structure increased and it looks like overall the smoothness is smoother and said than it is an X because these values are larger another useful way to actually look at the progress of an inversion is through a tea knopf curve and so what we're doing here is we plot Phi D against Phi M and so when we start the inversion we start with a very large Phi D and a very small fan because we started with our reference model and true model identical so if I am should be zero exactly but then as the inversion progresses we deke the value of Phi D well increasing the value of and so there's lots there's a number of criteria folks used to pick and choose what is your best inversion model the target misfit criteria would been P what's the first um what's the first inversion result below or target misfit other strategies try and look at where is the elbow of this kerf and we picked that model um and so you can you can play around with these so what I'm gonna do next is just define a couple more functions that are gonna allow us to sort of interactively explore our inversion results so in this case what we're doing is we're gonna plot there's three things we're plotting we've got the T knopf curve we've got our normalized misfit and we've got our resistivity model so at the initial iteration we really don't see much here our misfit basically looks like our data with plot because we haven't seen anything yet but as we sort of as we scroll through and we go from iteration to iteration you start to see a bit of structure introduced so we're seeing seeing this um conductor is starting to show up um and you can see here when we actually start to hit our target misfit so it looks like this literation perhaps we've we've gotten close to our target misfit so we can see both both blocks starting to show up now if we continue on the targets get stronger so perhaps if you wanted an inversion result more like this and if we were evaluating our noise model and things like that maybe we would re-evaluate our noise model if we think that this is is true in our target misfit was set to you high but if you keep going along what we can actually start to see is eventually there's gonna be some extra structure that starts to get introduced so we see some interesting features happening here and so what this is starting to show us is that our inversion is actually basically starting to go in and fits them in the noise and so even though we're pushing our data misfit really hard we're really really trying to respect the data that's not always something that you want to be doing because here we're introducing more structure than than what is actually geologically reasonable in this case and so that's why choosing choosing an appropriate as dropping criteria is an important thing to do okay so we're ten minutes out so I wanted to show you a couple more things that you can look at one thing that I did is I copied that little app code right here immediately below the inversion so that you can just run the inversion and explore note that this is sort of one of the features but it also can be a drawback if you forget about these things with the notebook is that this is run out of order right I need all of these functions in order for this to work so just remember if you comment this out or if you uncomment this and make it live and then you shut down this notebook and two weeks later you want to play around with it again you first need to run these cells below so there's a few different things you can go ahead and play with we can look at the Alpha values we can look at the influence of the noise floor and the reference model so the first ones I want to focus in on I'm just gonna reduce the maximum number of iterations so we can save ourselves a little bit of time um but let's say let's see what happens if I really increase alpha said let's say I know that this is I'm expecting vertical structures a pipe like structure and I just want alpha said to I want to really smooth things in this egg direction what happens so we can go ahead and rerun them and see what happens and we can explore explore those results so again we see we started at the same Phi D we started with the same model so that should be the case again Phi M is zero and now we can go ahead and look at these results and so in this case as we scroll through and look at our updates as a function of the iteration we can already see some very vertical structures being introduced this one's a little more subtle the resistors a bit more subtle but this conductor we're really seeing something come through quite vertical and if we look here this is our first iteration below the target misfit we're actually seeing some pretty strange structures I mean yes it does look like this is one coherent unit but we're starting to see a little bit of striping and so yeah so it's a slightly unrealistic geologic geologic structures but it is fitting Medina we flipped the data and this this model respects the data but for example if you were expecting some some smoothness in a given direction you can you can change these values may be a factor of two up to factor five ten um can still give you reasonable results so I encourage you to play around with these values and see sort of what structures can you get out that are geologically reasonable um by by playing around with these programs but but it's still giving you two just yet if we said these guys one other thing I want to show you because I think this is a useful useful understanding of basically how sensitive your data are is what happens if we change if we change the reference and the starting model so let's change the resistivity are starting in reference resistivity from the background model which we know is the true background iemon let's make it twice that choice is conductive so we've made this more conductive so we've reduced the resistivity by a factor of two so let's go ahead and run on version will again stay at ten iterations and we can see what happens so in this case when changing the the reference model like this this is a tool that's often used to estimate the depth of Investigation of an inversion result is running inversions with different starting and reference models and I'll show you why is because here if we then look at those results in our app so we'll scroll through so you can see that the colors different parents because we we started with a more conductive background and so as we go through iterations we can see we're starting to actually closer to the true background resistivity near the surface and we're starting to then see that resistive block and conductive block come through and then it looks like we're starting to fit the data around here here you could switch between those two inversion results fairly readily but we see so perhaps sort of like this sphere here a hemisphere data we're starting to get the right um background model but outside of that we're basically controlled by by what we started with and the regularization really is driving what we see there even if we continue sort of pushing the inversion as hard as we can these regions don't change and so what that tells you is that we don't really have any data that are telling us much about this region this region here is really controlled by the regularization and so our sensitivity the sensitivity of our data does not extend that far and so when folks talk about a depth of investigation this is one mechanism that you can use to try and understand which regions of the mesh and which regions of the domain do we actually have data that can support can support model changes and changes to the resistivity there and so I realized we are coming up on the hour and so what I want to leave you with with this inversion in this example is a few more ideas of things that you can share so as in the previous example in previous notebook what happens if you introduce a resistive layer a conductive layer how does that change our ability to resolve those blocks and then get creative let's see if you can build up an interesting model and try and stump some folks and try and try and see how hard you can see how well you can recover it and see what things are really challenging to recover in a DC resistivity inversion so please send us your interesting images and ideas on slack and we'd be happy to follow up and continue to conversation there and so with that I will say thank you to everyone they know pleasure being with you too and I hope you've learned some things I just wanted to put up the resources slide again so there's a few links in the GOC site for learning both about sort of the fundamentals of geophysical surveys as well as some aspects of inversion and then again to end with a shout-out to everyone I had the pleasure of speaking with you all today but this certainly would not have happened without an incredible team of people who've dedicated a lot of effort to simple and also specifically to to this tutorial so thanks to everyone here

Channel: Software Underground

Views: 7,825

Keywords: geophysics, inversion, modeling, programming, python, tutorial

Id: jZ7Sj9cnnso

Length: 180min 15sec (10815 seconds)

Published: Wed Jun 10 2020