In [1]:
from SimPEG import Mesh, EM, Utils, Maps
from matplotlib.colors import LogNorm
%pylab inline
import numpy as np
from scipy.constants import mu_0
from ipywidgets import interact, IntSlider
import cPickle as pickle
In [2]:
url = "https://storage.googleapis.com/simpeg/kevitsa_synthetic/"
files = ['dc_mesh.txt', 'dc_sigma.txt']
keys = ['mesh', 'sigma']
downloads = Utils.download([url + f for f in files], folder='./KevitsaDC', overwrite=True)
downloads = dict(zip(keys, downloads))
This model is a synthetic based on geologic surfaces interpreted from seismic data over the Kevitsa deposit in Finland. Synthetic 3D conductivity model is generated, and below figure shows conductivity section acrosses the mineralzined zone of interest. Nearsurface conductor on the lefthand side corresponds to sedimentary unit, and embedded conductor on the righthand side indicates conductive mineralized zone. Our interest here is in conductive mineralized zone at depth.
In [3]:
mesh3D = Mesh.TensorMesh.readUBC(downloads["mesh"])
sigmadc = mesh3D.readModelUBC(downloads["sigma"])
actind = ~np.isnan(sigmadc)
figsize(8, 4)
indy = 6
temp = 1./sigmadc.copy()
temp[~actind] = np.nan
out = mesh3D.plotSlice(temp, normal="Y", ind=indy, pcolorOpts={"norm": LogNorm(), "cmap":"jet_r"}, clim=(1e0, 1e3))
plt.ylim(-800, 250)
plt.xlim(5000, 11000)
plt.gca().set_aspect(2.)
plt.title(("y= %d m")%(mesh3D.vectorCCy[indy]))
cb = plt.colorbar(out[0], orientation="horizontal")
cb.set_label("Resistivity (Ohm-m)")
In [4]:
sig_halfspace = 2e-3
sig_target = 0.1
sig_air = 1e-8
times = np.logspace(-4, -2, 21)
def diffusion_distance(sigma, time):
return 1.28*np.sqrt(time/(sigma * mu_0))
print(
'min diffusion distance: {:.2e} m'.format(diffusion_distance(sig_halfspace, times.min()))
)
print(
'max diffusion distance: {:.2e} m'.format(diffusion_distance(sig_halfspace, times.max()))
)
# x-direction
csx = 20 # core mesh cell width in the x-direction
ncx = 20
npadx = 15 # number of x padding cells
# z-direction
csz = 20 # core mesh cell width in the z-direction
ncz = 40
npadz = 15 # number of z padding cells
# padding factor (expand cells to infinity)
pf = 1.3
# cell spacings in the x and z directions
hx = Utils.meshTensor([(csx, ncx), (csx, npadx, pf)])
hz = Utils.meshTensor([(csz, npadz, -pf), (csz, ncz), (csz, npadz, pf)])
# define a SimPEG mesh
mesh = Mesh.CylMesh([hx, 1, hz], x0 ="00C")
In [5]:
# X and Z limits we want to plot to. Try
xlim = np.r_[0., mesh.vectorCCx.max()]
zlim = np.r_[mesh.vectorCCz.max(), mesh.vectorCCz.min()]
fig, ax = plt.subplots(1,1)
mesh.plotGrid(ax=ax)
ax.set_title('Simulation Mesh')
ax.set_xlim(xlim)
ax.set_ylim(zlim)
print(
'The maximum diffusion distance (in background) is: {:.2e} m. '
'Does the mesh go sufficiently past that?'.format(
diffusion_distance(sig_halfspace, times.max())
)
)
ax.set_aspect("equal")
Next, we put the model on the mesh
In [6]:
# create a vector that has one entry for every cell center
sigma = sig_air*np.ones(mesh.nC) # start by defining the conductivity of the air everwhere
sigma[mesh.gridCC[:,2] < 0.] = sig_halfspace # assign halfspace cells below the earth
sigma_background = sigma.copy()
sigma_layer = sigma.copy()
radius = 150.
# indices of the sphere (where (x-x0)**2 + (z-z0)**2 <= R**2)
layer_ind = np.logical_and(mesh.gridCC[:,2]>-300, mesh.gridCC[:,2]<-200)
blk_ind = (mesh.gridCC[:,0] < radius) & layer_ind
sigma[blk_ind] = sig_target # assign the conductivity of the sphere
sigma_layer[layer_ind] = sig_target # assign the conductivity of the sphere
In [7]:
plt.set_cmap(plt.get_cmap('jet_r'))
# Plot a cross section of the conductivity model
fig, ax = plt.subplots(1,1)
out = mesh.plotImage(np.log10(1./sigma_layer), ax=ax, mirror=True, clim=(0, 3), grid=False)
cb = plt.colorbar(out[0], ticks=np.linspace(0,3,4), format="10$^{%.1f}$")
# plot formatting and titles
cb.set_label('Resistivity (Ohm-m)', fontsize=13)
ax.axis('equal')
ax.set_xlim([-120., 120.])
ax.set_ylim([-500., 0.])
ax.set_title('Layer')
Out[7]:
In [8]:
# Plot a cross section of the conductivity model
fig, ax = plt.subplots(1,1)
out = mesh.plotImage(np.log10(1./sigma), ax=ax, mirror=True, clim=(0, 3), grid=False)
# plot formatting and titles
cb = plt.colorbar(out[0], ticks=np.linspace(0,3,4), format="10$^{%.1f}$")
# plot formatting and titles
cb.set_label('Resistivity (Ohm-m)', fontsize=13)
ax.axis('equal')
ax.set_xlim([-120., 120.])
ax.set_ylim([-500., 0.])
ax.set_title('Cylinder')
Out[8]:
In [9]:
rx_loc = np.array([[0., 0., 41.]])
src_loc = np.array([[0., 0., 41.]])
offTime = 0.007307
peakTime = 0.006
a = 3.
dbdt_z = EM.TDEM.Rx.Point_dbdt(locs=rx_loc, times=times+offTime, orientation='z') # vertical db_dt
rxList = [dbdt_z] # list of receivers
srcList = [
EM.TDEM.Src.CircularLoop(
rxList, loc=src_loc, radius=13., orientation='z', waveform=EM.TDEM.Src.VTEMWaveform(offTime=offTime, peakTime=peakTime, a=3.)
)
]
In [10]:
# solve the problem at these times
timeSteps = [(peakTime/5, 5), ((offTime-peakTime)/5, 5), (1e-5, 10), (5e-5, 10), (1e-4, 10), (5e-4, 19)]
prob = EM.TDEM.Problem3D_b(mesh, timeSteps = timeSteps, sigmaMap=Maps.IdentityMap(mesh))
survey = EM.TDEM.Survey(srcList)
prob.pair(survey)
In [11]:
src = srcList[0]
rx = src.rxList[0]
wave = []
for time in prob.times:
wave.append(src.waveform.eval(time))
wave = np.hstack(wave)
In [12]:
plt.plot(prob.times, wave, 'k.-')
plt.plot(rx.times, np.zeros_like(rx.times), 'r.')
plt.ylim(-0.2, 1.2)
plt.grid(True)
plt.title('Current Waveform')
plt.xlabel('time (s)')
Out[12]:
In [13]:
d_background = survey.dpred(sigma_background)
d_layer = survey.dpred(sigma_layer)
d = survey.dpred(sigma)
In [14]:
area = 13**2*np.pi
figsize(6, 3)
plt.loglog((rx.times-offTime)*1e6, -d_layer*1e12/area, 'k', lw=2)
plt.loglog((rx.times-offTime)*1e6, -d*1e12/area , 'b', lw=2)
plt.loglog((rx.times-offTime)*1e6, -d_background*1e12/area, 'k--', lw=1)
plt.xlabel("Time (micro-s)")
plt.ylabel("Voltage (pV/A-m$^4$)")
plt.legend(("Layer", "Cylinder","Half-space"), loc=1, fontsize = 10)
plt.ylim(1e-4, 1e1)
plt.grid(True)
What was your thoughts on the above plot? Can we see conductive mineralzied zone? Singals from Layer and Cylinder have significant difference, can you explain why?
Underlying physics of the measured voltage can be govered by Faraday's law:
$$ \nabla \times \vec{e} - \frac{d\vec{b}}{dt}$$By showing how electric field propagates in the subsurface we illustrate why layer and cylinder model show significant difference.
In [15]:
f_layer = prob.fields(sigma_layer)
In [16]:
plt.set_cmap(plt.get_cmap('viridis'))
def vizfield_layer(itime):
fig = plt.figure(figsize = (7*0.8,5*0.8))
ax = plt.subplot(111)
cb = plt.colorbar(mesh.plotImage(mesh.aveE2CC*f_layer[src, 'e', itime], ax=ax, mirror=True)[0])
# plot formatting and titles
cb.set_label('e$_{y}$ (V/m)', fontsize=13)
ax.axis('equal')
ax.set_xlim([-300., 300.])
ax.set_ylim([-500., 0.])
ax.set_title(('|e$_{y}$| at %d micro-s')%(prob.times[itime]*1e6))
plt.show()
interact(vizfield_layer, itime=IntSlider(min=0, max=len(prob.times)-1, step=1, value=11))
Out[16]:
In [17]:
f = prob.fields(sigma)
In [18]:
def vizfield_cylinder(itime):
fig = plt.figure(figsize = (7*0.8,5*0.8))
ax = plt.subplot(111)
cb = plt.colorbar(mesh.plotImage(mesh.aveE2CC*f[src, 'e', itime], ax=ax, mirror=True)[0])
# plot formatting and titles
cb.set_label('e$_{y}$ (V/m)', fontsize=13)
# ax.axis('equal')
ax.set_xlim([-300., 300.])
ax.set_ylim([-500., 0.])
ax.set_title(('|e$_{y}$| at %d micro-s')%(prob.times[itime]*1e6))
plt.tight_layout()
plt.show()
interact(vizfield_cylinder, itime=IntSlider(min=0, max=len(prob.times)-1, step=1, value=11))
Out[18]: