In [27]:
from SimPEG import *
%matplotlib inline
import matplotlib.pyplot as plt
DC resistivity is a zero-frequency electromagnetic problem. A DC resistivity geophysical survey is sensitive to variations in electrical conductivity. To conduct a DC resistivity survey, we use a source source electrode to inject a current into the subsurface and sink to pick up that current.
The governing equation is:
$$ \nabla \cdot \sigma \nabla \phi = -s $$where $\sigma$ is the electrical conductivity (S/m) which varies depending on the subsurface materials, $\phi$ is the electric potential (V), and $s$ is the source term, which captures the impact of both the source and sink.
We will solve the DC problem using a finite volume simulation, which requires that we set up a mesh on which to compute. For this example, we use a 2D mesh, but this can be expanded to 3D using mesh = Mesh.TensorMesh([nx,ny,nz]).
We use the SimPEG package (http://simpeg.xyz/) to perform the simulation.
In [28]:
# Define a unit square domain
# This can be hidden and imported depending on how much you want to show
nx,ny = 60,60 # number of cells in x,y
mesh = Mesh.TensorMesh([nx,ny]) # build a tensor mesh
sigma = np.ones(mesh.nC) # assign a conductivity model
# create source
xp, yp = 0.25, 0.5
xn, yn = 0.75, 0.5
sigmaback = 1.
s = sigmaback*np.zeros(mesh.nC)
indp = Utils.closestPoints(mesh,np.r_[xp,yp],'CC')
indn = Utils.closestPoints(mesh,np.r_[xn,yn],'CC')
s[indp] = 1.
s[indn] = -1
# conductive block
sigmabl = 100.
p0 = np.r_[0.4,0.4]
p1 = np.r_[0.6,0.6]
In [29]:
# Build differential operators using the SimPEG mesh class
# This can be hidden and imported depending on how much you want to expose
def getOperators(mesh,sigma):
Div = mesh.faceDiv
Sigma = Utils.sdiag(1/(mesh.aveF2CC.T * (1/sigma)))
Grad = mesh.cellGrad
return Div, Sigma, Grad
def getA(mesh,sigma):
Div, Sigma, Grad = getOperators(mesh,sigma)
return Div * Sigma * Grad
In [30]:
# Construct A Matrix
Div, Sigma, Grad = getOperators(mesh,sigma)
A = Div * Sigma * Grad # looks like the equation!
Ainv = Solver(A)
In [31]:
phi = Ainv * -s # solve
In [32]:
mesh.plotImage(phi)
plt.title('Electric Potential, $\phi$')
Out[32]:
In [33]:
# this can be hidden and imported
def SolveDC(xp,yp,xn,yn,mesh,sigma):
# set up the source
s = np.zeros(mesh.nC)
indp = Utils.closestPoints(mesh,np.r_[xp,yp],'CC')
indn = Utils.closestPoints(mesh,np.r_[xn,yn],'CC')
s[indp] = 1.
s[indn] = -1.
# get A
A = getA(mesh,sigma)
A[0,0] *= 1/mesh.vol[0] # correct for the null space
Ainv = Solver(A)
# solve
phi = Ainv * -s
phi -= phi[-1]
return phi
def plotPhi(phi):
mesh.plotImage(phi)
plt.title('Electric Potential, $\phi$')
plt.show()
The model has a uniform background conductivity:
sigma = 1 S/mWe excite the system using a source which has positive and negative electrodes (ie. a current source and a current). We will plot the electric potential; when a DC resistivity survey is conducted, the data are potential differences, measured in units of volts.
The parameters you can explore here are:
xp: x-position of the positive electrodeyp: y-position of the positive electrodexn: x-position of the negative electrodeyn: y-position of the negative electrode
In [34]:
from IPython.html.widgets import interact
interact(lambda xp,yp,xn,yn: plotPhi(SolveDC(xp,yp,xn,yn,mesh,sigma)), xp=[0.1,0.8,0.05],yp=[0.1,0.9,0.05],xn=[0.2,0.9,0.05],yn=[0.1,0.9,0.05])
In [35]:
def CondBlock(xp,yp,xn,yn,mesh,sigma,sigmabl,p0,p1):
sigmaBl = Utils.ModelBuilder.defineBlock(mesh.gridCC, p0, p1, np.r_[sigmabl,sigmaback])
return SolveDC(xp,yp,xn,yn,mesh,sigmaBl)
What happens if we add a conductive block?
A conductive block with a conductivity of 100 S/m in the 1 S/m background. If you move the positive and negative electrodes around, can you find it? How big is it?
In [36]:
interact(lambda xp,yp,xn,yn: plotPhi(CondBlock(xp,yp,xn,yn,mesh,sigma,sigmabl,p0,p1)), xp=[0.1,0.8,0.05],yp=[0.1,0.9,0.05],xn=[0.2,0.9,0.05],yn=[0.1,0.9,0.05])
Out[36]:
Write a function to compute the secondary field
In [37]:
def CondBlockSecondary(xp,yp,xn,yn,mesh,sigma,sigmabl,p0,p1):
sigmaBl = Utils.ModelBuilder.defineBlock(mesh.gridCC, p0, p1, np.r_[sigmabl,sigmaback])
sec = SolveDC(xp,yp,xn,yn,mesh,sigmaBl) - SolveDC(xp,yp,xn,yn,mesh,sigma)
return sec
Here, we plot the difference in potentials computed with and without the block, ie. the signal due to the presence of the block
In [38]:
interact(lambda xp,yp,xn,yn: plotPhi(CondBlockSecondary(xp,yp,xn,yn,mesh,sigma,sigmabl,p0,p1)), xp=[0.1,0.8,0.05],yp=[0.1,0.9,0.05],xn=[0.2,0.9,0.05],yn=[0.1,0.9,0.05])
Here, you can play with all of the parameters in this model:
xp: x-position of the positive electrodeyp: y-position of the positive electrodexn: x-position of the negative electrodeyn: y-position of the negative electrodelog10_sigmabl: $log_{10}(\sigma_{\text{block}})$blcx : x-center of the blockblcy : y-center of the blockbldx : x-width of the blockbldy : y-width of the blockThe first widget shows the measured potential and the second shows the secondary potential
In [39]:
# total potential
interact(lambda xp,yp,xn,yn,log10_sigmabl,blcx,blcy,bldx,bldy: plotPhi(CondBlock(xp,yp,xn,yn,mesh,sigma,10.**(log10_sigmabl),np.r_[blcx-bldx/2.,blcy-bldy/2.],np.r_[blcx+bldx/2.,blcy+bldy/2.])), xp=[0.1,0.8,0.05],yp=[0.1,0.9,0.05],xn=[0.2,0.9,0.05],yn=[0.1,0.9,0.05],log10_sigmabl=[-5,5,0.5],blcx=[0.,1.,0.05],blcy=[0.,1.,0.05],bldx=[0.,1.,0.05],bldy=[0.,1.,0.05])
Out[39]:
In [40]:
# secondary potential
interact(lambda xp,yp,xn,yn,log10_sigmabl,blcx,blcy,bldx,bldy: plotPhi(CondBlockSecondary(xp,yp,xn,yn,mesh,sigma,10.**(log10_sigmabl),np.r_[blcx-bldx/2.,blcy-bldy/2.],np.r_[blcx+bldx/2.,blcy+bldy/2.])), xp=[0.1,0.8,0.05],yp=[0.1,0.9,0.05],xn=[0.2,0.9,0.05],yn=[0.1,0.9,0.05],log10_sigmabl=[-5,5,0.5],blcx=[0.,1.,0.05],blcy=[0.,1.,0.05],bldx=[0.,1.,0.05],bldy=[0.,1.,0.05])
Out[40]:
In [ ]:
In [ ]: