Keywords: DC resistivity, forward simulation, 2.5D, apparent resistivity, tree mesh.
Summary: Here we use the module simpeg
Learning Objectives:
- How to define DC resistivity lines manually or by using utility functions in SimPEG.
- How to design a 2D tree mesh (or tensor mesh) for accurately simulating DC resistivity data for the 2.5D approach.
- How to define the Earth’s electrical properties according to conductivity OR resistivity.
- How to included surface topography in the forward simulation.
- How to plot simulated data in pseudosection.
Import Modules¶
Here, we import all of the functionality required to run the notebook for the tutorial exercise. All of the functionality specific to DC resistivity is imported from simpeg
# SimPEG functionality
from simpeg.electromagnetics.static import resistivity as dc
from simpeg.utils import model_builder
from simpeg.utils.io_utils.io_utils_electromagnetics import write_dcip2d_ubc
from simpeg import maps, data
from simpeg.electromagnetics.static.utils.static_utils import (
generate_dcip_sources_line,
pseudo_locations,
plot_pseudosection,
apparent_resistivity_from_voltage,
)
# discretize functionality
from discretize import TreeMesh
from discretize.utils import active_from_xyz
# Common Python functionality
import os
import numpy as np
from scipy.interpolate import interp1d
import matplotlib as mpl
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
mpl.rcParams.update({"font.size": 14}) # default font size
write_output = False # Optional
Define the Topography¶
True surface topography is defined as an (N, 3) numpy.ndarray. For use in a 2.5D simulation however, topography is defined as an (N, 2) numpy.ndarray, where the first coordinate represent along-line position and the second coordinate represents the vertical position. Here, we define the 2D topography used for the 2.5D simulation.
# Along-line locations
x_topo = np.linspace(-2000, 2000, 401)
# Elevation as a function of along-line location
T = 800.0
z_topo = 20.0 * np.sin(2 * np.pi * x_topo / T) + 140.0
z_topo[x_topo < -3 * T / 4] = 160.0
z_topo[x_topo > 3 * T / 4] = 120.0
z_topo += 50.0 * (1.0 + np.tanh(-3 * (x_topo + 1200.0) / T))
z_topo -= 50.0 * (1.0 + np.tanh(3 * (x_topo - 1200.0) / T))
# Define full 2D topography
topo_2d = np.c_[x_topo, z_topo]
# Plot 2D topography
fig = plt.figure(figsize=(10, 2))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
ax.plot(x_topo, z_topo, color="b", linewidth=2)
ax.set_xlabel("x (m)", labelpad=5)
ax.set_ylabel("z (m)", labelpad=5)
ax.grid(True)
ax.set_title("Topography (Exaggerated z-axis)", fontsize=16, pad=10)
plt.show(fig)
Define the Survey¶
DC (and IP) surveys within SimPEG require the user to create and connect three types of objects:
- receivers: which defines the locations of the potential (or MN) electrodes and the type of data; e.g. ‘volt’ for normalized voltage (V/A), ‘apparent_resistivity’ for apparent resistivity () or ‘apparent_chargeability’ for apparent chargeability (unitless). Note only M electrode locations are needed to define pole receivers.
- sources: which defines the locations of the current (or AB) electrodes, and their associated receivers. Note only A electrode locations are needed to define pole sources.
- survey: the object which stores and organizes all of the sources and receivers.
Each DC/IP datum corresponds to a unique pair of current and potential electrode locations. Within the SimPEG framework, there is a standard way in which the electrodes are organized. Because each current source (pole or dipole) is discretized to form a separate right-hand side for the DC resistivity PDE, each source object defines the electrodes for a single source. But for each source object, we can instantiate a single receiver object to store the electrode locations of all potential electrodes.
Using SimPEG, there are three approaches one might use to generate the a 2.5D resistivity survey. The final approach will be used to simulate the data for this tutorial.
Option A: Define each source and its associated receivers directly
For a 2.5D simulation, current electrode locations are defined as (2,) numpy.array and the associated set of potential electrode locations are defined as (N, 2) numpy.ndarray. We can define Pole or Dipole sources. And we can define Pole or Dipole receivers.
Here, the survey consists of an 800 m long EW dipole-dipole line with an electrode spacing of 40 m. There is a maximum of 10 potential electrodes per current electrode. And the data defined by the receivers are current-normalized voltages in V/A.
# Define survey line parameters
survey_type = "dipole-dipole"
dimension_type = "2D"
data_type = "volt"
end_locations = np.r_[-400.0, 400.0] # along-line position
station_separation = 40.0
num_rx_per_src = 10
# Define linear interpolation function for elevation
interp_fun = interp1d(x_topo, z_topo)
# Define electrode locations
electrode_locations_x = np.arange(
end_locations[0], end_locations[1] + station_separation, station_separation
)
electrode_locations_z = interp_fun(electrode_locations_x)
electrode_locations = np.c_[electrode_locations_x, electrode_locations_z]
# Number of electrode locations
n_electrodes = len(electrode_locations_x)
# Instantiate empty list for sources
source_list = []
ii = 0
while ii < n_electrodes - 3:
# A and B electrode locations
location_a = electrode_locations[ii, :]
location_b = electrode_locations[ii + 1, :]
# M and N electrode locations
ii_max = np.min([ii + 3 + num_rx_per_src, n_electrodes])
locations_m = electrode_locations[ii + 2 : ii_max - 1]
locations_n = electrode_locations[ii + 3 : ii_max]
# Define receivers for source ii
receivers_list = [
dc.receivers.Dipole(
locations_m=locations_m, locations_n=locations_n, data_type=data_type
)
]
# Append source ii to list
source_list.append(
dc.sources.Dipole(receivers_list, location_a=location_a, location_b=location_b)
)
ii += 1
# Define survey
survey = dc.Survey(source_list)
Option B: Survey from ABMN electrode locations
If we have (N, 2) numpy.ndarray for A, B, M and N electrode locations for each datum (loaded or created), we can use the generate
Option C: Survey from a set of survey lines
If the survey is comprised of a single DC resistivity line, we can use the generate
# Generate source list for DC survey line
source_list = generate_dcip_sources_line(
survey_type,
data_type,
dimension_type,
end_locations,
topo_2d,
num_rx_per_src,
station_separation,
)
# Define survey
survey = dc.survey.Survey(source_list, survey_type=survey_type)
We can extract various objects and properties from the objects used to generate the survey. As we can see, all receivers associated with each source are defined within a single object.
fig = plt.figure(figsize=(10, 2))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])
ax.plot(x_topo, z_topo, color="b", linewidth=1)
ax.scatter(electrode_locations_x, electrode_locations_z, 8, "r")
ax.set_xlim([x_topo.min(), x_topo.max()])
ax.set_xlabel("x (m)", labelpad=5)
ax.set_ylabel("z (m)", labelpad=5)
ax.grid(True)
ax.set_title("Topography and electrode locations", fontsize=16, pad=10)
plt.show(fig)
Here we plot the electrode locations. We use the pseudo_locations utility to extract the pseudo-locations.
pseudo_locations_xz = pseudo_locations(survey)
fig = plt.figure(figsize=(8, 2.75))
ax = fig.add_axes([0.1, 0.1, 0.85, 0.8])
ax.scatter(pseudo_locations_xz[:, 0], pseudo_locations_xz[:, -1], 8, "r")
ax.set_xlabel("x (m)")
ax.set_ylabel("z (m)")
ax.set_title("Pseudo-locations")
plt.show()
Design a (Tree) Mesh¶
Meshes are designed using the discretize package. See the discretize user tutorials to learn more about creating meshes. Here, the forward simulation is computed for a tree mesh. Because of the modular nature of SimPEG, you could define a tensor mesh instead.
Standard approach for DC/IP: The electric potential produced by a current electrode falls off as . So smaller cells are needed near the current electrodes to model the fields accurately, and larger cells can be used away from the current electrodes where the fields are smooth. Tree meshes are well-suited for DC (and IP) simulation because the cell size can be increased at specified distances from the current electrodes. For DC resistivity meshing, we advise the following considerations and rules of thumb:
- Because there are no currents in the air, we do not need to pad upwards. I.e. the top of the mesh corresponds to the top of the topography.
- We require at least 2-3 cells between each current electrode; with more accurate results being obtained when the minimum cell size is smaller. For a 2.5D problem geometry, we can discretize much finer.
- To be safe, the padding thickness should be at least 2-3 times the largest electrode spacing.
- The increase in cell size at increasing distances from the current electrodes should not happen too abruptly. At each cell size, you should have a layer at least 4 cells thick before increasing the cell size.
- Finer discretization is required when topography is significant.
Tutorial mesh: Here, a minimum cell width of 4 m (or 1/10 the minimum electrode spacing) is used within our survey region. The largest electrode spacing was 400 m, so a the padding was extended at least 1200 m from the survey region. Using the refine_surface method, we refine the tree mesh where there is significant topography. And using the refine_points methods, we refine based on electrodes locations. Visit the tree mesh API to see additional refinement methods.
dh = 4 # base cell width
dom_width_x = 3200.0 # domain width x
dom_width_z = 2400.0 # domain width z
nbcx = 2 ** int(np.round(np.log(dom_width_x / dh) / np.log(2.0))) # num. base cells x
nbcz = 2 ** int(np.round(np.log(dom_width_z / dh) / np.log(2.0))) # num. base cells z
# Define the base mesh with top at z = 0 m.
hx = [(dh, nbcx)]
hz = [(dh, nbcz)]
mesh = TreeMesh([hx, hz], x0="CN", diagonal_balance=True)
# Shift top to maximum topography
mesh.origin = mesh.origin + np.r_[0.0, z_topo.max()]
# Mesh refinement based on topography
mesh.refine_surface(
topo_2d,
padding_cells_by_level=[0, 0, 4, 4],
finalize=False,
)
# Extract unique electrode locations.
unique_locations = survey.unique_electrode_locations
# Mesh refinement near electrodes.
mesh.refine_points(
unique_locations, padding_cells_by_level=[8, 12, 6, 6], finalize=False
)
mesh.finalize()
If desired, we can extract various properties of the mesh. E.g.
print("# of cells: {}".format(mesh.n_cells)) # Number of cells
print("# of x-faces: {}".format(mesh.n_faces_x)) # Number of x-faces
print("Origin: {}".format(mesh.origin)) # bottom-southewest corner
print("Max cell volume: {}".format(mesh.cell_volumes.max())) # Largest cell size
print(mesh.cell_centers[0:5, :]) # Cell center locations
# of cells: 12311
# of x-faces: 12224
Origin: [-2048. -1788.24726232]
Max cell volume: 262144.0
[[-1792. -1532.24726232]
[-1280. -1532.24726232]
[-1792. -1020.24726232]
[-1280. -1020.24726232]
[ -768. -1532.24726232]]
And we can plot the mesh and electrode locations.
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_axes([0.14, 0.17, 0.8, 0.7])
mesh.plot_grid(ax=ax1, linewidth=1)
ax1.grid(False)
ax1.set_xlim(-1500, 1500)
ax1.set_ylim(np.max(z_topo) - 1000, np.max(z_topo))
ax1.set_title("Mesh")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("z (m)")
plt.show()
Define the Active Cells¶
Simulated geophysical data are dependent on the subsurface distribution of physical property values. As a result, the cells lying below the surface topography are commonly referred to as ‘active cells’. And air cells, whose physical property values are fixed, are commonly referred to as ‘inactive cells’. Here, the discretize active_from_xyz utility function is used to find the indices of the active cells using the mesh and surface topography. The output quantity is a bool
array.
# Indices of the active mesh cells from topography (e.g. cells below surface)
active_cells = active_from_xyz(mesh, topo_2d)
# number of active cells
n_active = np.sum(active_cells)
Models and Mappings¶
In SimPEG, the term ‘model’ is not necessarily synonymous with the set of physical property values defined on the mesh. For example, the model may be defined as the logarithms of the physical property values, or be the parameters defining a layered Earth geometry. Models in SimPEG are 1D numpy.ndarray whose lengths are equal to the number of model parameters.
Classes within the simpeg.maps module are used to define the mapping that connects the model to the physical property values used in the DC resistivity simulation. Sophisticated mappings can be defined by combining multiple mapping objects. But in the simplest case, the mapping is an identity map and the model consists of the conductivity/resistivity values for all mesh cells (including air).
When simulating DC resistivity data, we have the choice of using resistivity or conductivity to define the Earth’s electrical properties. Here, we define the model and its associate mapping for two cases:
- The model consists of the conductivity values for all active cells
- The model consists of the log-resistivity values for all active cells
Define the Model¶
The units for resistivity are and the units for conductivity are .
air_conductivity = 1e-8
background_conductivity = 1e-2
conductor_conductivity = 1e-1
resistor_conductivity = 1e-3
# Define conductivity model
conductivity_model = background_conductivity * np.ones(n_active)
ind_conductor = model_builder.get_indices_sphere(
np.r_[-120.0, 40.0], 60.0, mesh.cell_centers[active_cells, :]
)
conductivity_model[ind_conductor] = conductor_conductivity
ind_resistor = model_builder.get_indices_sphere(
np.r_[120.0, 72.0], 60.0, mesh.cell_centers[active_cells, :]
)
conductivity_model[ind_resistor] = resistor_conductivity
# Define log-resistivity model
log_resistivity_model = np.log(1 / conductivity_model)
Mapping from the Model to the Mesh¶
For our first case, we use the simpeg
For the second case, we both the simpeg
# Conductivity map. Model parameters are conductivities for all active cells.
conductivity_map = maps.InjectActiveCells(mesh, active_cells, air_conductivity)
# Resistivity map. Model parameters are log-resistivities for all active cells.
log_resistivity_map = maps.InjectActiveCells(
mesh, active_cells, 1 / air_conductivity
) * maps.ExpMap(nP=n_active)
Plot the Model¶
To show the geometry of the problem, we plot the conductivity model using the plot_slice method.
# Generate a mapping to ignore inactice cells in plot
plotting_map = maps.InjectActiveCells(mesh, active_cells, np.nan)
fig = plt.figure(figsize=(9, 4))
norm = LogNorm(vmin=1e-3, vmax=1e-1)
ax1 = fig.add_axes([0.14, 0.17, 0.68, 0.7])
mesh.plot_image(
plotting_map * conductivity_model,
ax=ax1,
grid=False,
pcolor_opts={"norm": norm, "cmap": mpl.cm.RdYlBu_r},
)
ax1.set_xlim(-500, 500)
ax1.set_ylim(-300, 200)
ax1.set_title("Conductivity Model")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("z (m)")
ax2 = fig.add_axes([0.84, 0.17, 0.03, 0.7])
cbar = mpl.colorbar.ColorbarBase(
ax2, norm=norm, orientation="vertical", cmap=mpl.cm.RdYlBu_r
)
cbar.set_label(r"$\sigma$ (S/m)", rotation=270, labelpad=15, size=12)
plt.show()
Project Electrodes to Discretized Topography¶
Surface DC resistivity data will not be modeled accurately if the electrodes are modeled as living above or below the surface. It is especially problematic when electrodes are modeled as living in the air. Prior to simulating surface DC resistivity data, we must project the electrodes from their true elevation to the surface of the discretized topography. This is done using the drape
survey.drape_electrodes_on_topography(mesh, active_cells, option="top")
Define the Forward Simulation¶
In SimPEG, the physics of the forward simulation is defined by creating an instance of an appropriate simulation class. There are two simulation classes which may be used to simulate 2.5D DC resistivity data:
- Simulation2DNodel, which defines the discrete electric potentials on mesh nodes.
- Simulation2DCellCentered, which defines the discrete electric potentials at cell centers.
For surface DC resistivity data, the nodal formulation is more well-suited and will be used here. The cell-centered formulation works well for simulating borehole DC resistivity data. To fully define the forward simulation, we need to connect the simulation object to:
- the survey
- the mesh
- the mapping from the model to the mesh
This is accomplished by setting each one of the aforementioned items as a property of the simulation object. Here, we define two simulation objects, one where the model defines the subsurface conductivities, and one where the model defines subsurface log-resistivities. When our model is used to define subsurface electric conductivity, the mapping is set using the sigmaMap
keyword argument. However when our model is used to define subsurface electric resistivity, the mapping must be set using the rhoMap
keyword argument
# DC simulation for a conductivity model
simulation_con = dc.simulation_2d.Simulation2DNodal(
mesh, survey=survey, sigmaMap=conductivity_map
)
# DC simulation for a log-resistivity model
simulation_res = dc.simulation_2d.Simulation2DNodal(
mesh, survey=survey, rhoMap=log_resistivity_map
)
Predict DC Resistivity Data¶
Once any simulation within SimPEG has been properly constructed, simulated data for a given model vector can be computed using the dpred method. Note that despite the difference in how we defined the models representing the Earth’s electrical properties, the data predicted by both simulations is equivalent.
dpred_con = simulation_con.dpred(conductivity_model)
/t40array/ssoler/miniforge3/envs/simpeg-user-tutorials/lib/python3.10/site-packages/simpeg/electromagnetics/static/resistivity/simulation_2d.py:756: RuntimeWarning: invalid value encountered in divide
r_hat = r_vec / r[:, None]
/t40array/ssoler/miniforge3/envs/simpeg-user-tutorials/lib/python3.10/site-packages/simpeg/electromagnetics/static/resistivity/simulation_2d.py:783: RuntimeWarning: invalid value encountered in divide
alpha[not_top] = (ky * k1e(ky * r) / k0e(ky * r) * r_dot_n)[not_top]
/t40array/ssoler/miniforge3/envs/simpeg-user-tutorials/lib/python3.10/site-packages/simpeg/simulation.py:197: DefaultSolverWarning: Using the default solver: Pardiso.
If you would like to suppress this notification, add
warnings.filterwarnings('ignore', simpeg.utils.solver_utils.DefaultSolverWarning)
to your script.
return get_default_solver(warn=True)
dpred_res = simulation_res.dpred(log_resistivity_model)
print("MAX ABSOLUTE ERROR = {}".format(np.max(np.abs(dpred_con - dpred_res))))
MAX ABSOLUTE ERROR = 1.970645868709653e-15
Plot Data in Pseudosection¶
Here we use the plot_pseudosection utility function to represent the predicted data on a pseudosection plot as apparent conductivities. Since our receivers were defined to simulate data as normalized voltages, we use the apparent
# Plot voltages pseudo-section
fig = plt.figure(figsize=(8, 2.75))
ax1 = fig.add_axes([0.1, 0.15, 0.75, 0.78])
plot_pseudosection(
survey,
dobs=np.abs(dpred_con),
plot_type="scatter",
ax=ax1,
scale="log",
cbar_label="V/A",
scatter_opts={"cmap": mpl.cm.viridis},
)
ax1.set_title("Normalized Voltages")
plt.show()
# Get apparent conductivities from volts and survey geometry
apparent_conductivities = 1 / apparent_resistivity_from_voltage(survey, dpred_con)
# Plot apparent conductivity pseudo-section
fig = plt.figure(figsize=(8, 2.75))
ax1 = fig.add_axes([0.1, 0.15, 0.75, 0.78])
plot_pseudosection(
survey,
dobs=apparent_conductivities,
plot_type="contourf",
ax=ax1,
scale="log",
cbar_label="S/m",
mask_topography=True,
contourf_opts={"levels": 20, "cmap": mpl.cm.RdYlBu_r},
)
ax1.set_title("Apparent Conductivity")
plt.show()
Optional: Write DC resistivity data and topography.
if write_output:
dir_path = os.path.sep.join([".", "fwd_dcr_2d_outputs"]) + os.path.sep
if not os.path.exists(dir_path):
os.mkdir(dir_path)
# Add 5% Gaussian noise to each datum
rng = np.random.default_rng(seed=225)
std = 0.05 * np.abs(dpred_con)
dc_noise = rng.normal(scale=std, size=len(dpred_con))
dobs = dpred_con + dc_noise
# Create a survey with the original electrode locations
# and not the shifted ones
# Generate source list for DC survey line
source_list = generate_dcip_sources_line(
survey_type,
data_type,
dimension_type,
end_locations,
topo_2d,
num_rx_per_src,
station_separation,
)
survey_original = dc.survey.Survey(source_list)
# Write out data at their original electrode locations (not shifted)
data_obj = data.Data(survey_original, dobs=dobs, standard_deviation=std)
fname = dir_path + "dc_data.obs"
write_dcip2d_ubc(fname, data_obj, "volt", "dobs")
fname = dir_path + "topo_2d.txt"
np.savetxt(fname, topo_2d, fmt="%.4e")