ERTLab is the 3D resistivity and chargeability inversion software that has. For the design of the electrode geometry, resistivity/chargeability forward modeling,. Resistivity and IP data quality control, denoising and editing data module, EM removing. A set of tools for mathematical modeling of geoelectrical section. Standard mesh “bitmap” style, arbitrary layered model and polygonal style model. Additional tools for gravity.
AbstractA Python package, pyres, was written to handle common as well as specialized input and output tasks for the R2 electrical resistivity (ER) modeling program. Input steps including handling field data, creating quadrilateral or triangular meshes, and data filtering allow repeatable and flexible ER modeling within a programming environment. Pyres includes non-trivial routines and functions for locating and constraining specific known or separately-parameterized regions in both quadrilateral and triangular meshes. Three basic examples of how to run forward and inverse models with pyres are provided. The importance of testing mesh convergence and model sensitivity are also addressed with higher-level examples that show how pyres can facilitate future research-grade ER analyses.
The basic workflow for pyres follows the standard workflow for running R2 without Python. To run a forward model in R2, the finite element mesh must be defined with either node spacing information or from a pre-made mesh file, along with the electrode combinations and their locations, an earth model that provides the ER structure that is to be investigated, and numerous model options built into R2.
For an inverse model, additional model options are specified for the inverse problem along with providing the ER measurements for each electrode combination, in addition to the inputs required for the forward model.Given these steps for modeling ER data in R2, pyres first prescribes the electrode configuration either by defining the number and spacing of electrodes or by loading in survey data through pyresutils.loaderdata. Next, the discretization of the mesh is defined along with the definition of the length scales controlling the size of the foreground and background solution domains. If the mesh is made of quadrilateral elements, then the discretization of the mesh nodes can be defined using NumPy arrays. For a triangular mesh, pyres includes functionality to construct a triangular finite element mesh using the freeware program Gmsh using the meshtools.meshR2 class. Regions within the mesh with a priori information are defined either during the mesh creation for a triangular mesh or after the mesh arrays are made for a quadrilateral mesh (see section.
R2 examples). Since pyres constructs the mesh, the node and element indices required for assigning regions of a priori information or constrained with additional parameterization can be readily located and incorporated into the R2 input data files. Similarly, topographic information can be added to quadrilateral and triangular meshes (see R2 example 2 for quadrilateral mesh). Next, the two input files for R2, protocol.dat and R2.in, are created with the r2 tools.R2 class using an electrode configuration defined by measured input data or a synthetic survey with the pyresutils.makesyntheticsurvey function. The r2 tools.R2 class can then be used to directly run R2 with the r2 tools.R2.runr2 command without copying the executable into the directory containing the R2 input files.In addition to managing R2 inputs and running R2 forward and inverse models, pyres also contains several other input and output methods that can be incorporated with user-defined Python code for higher-level pre- or post-processing steps. When working with noisy ER data, pseudo-sections plotting apparent resistivity as a function of dipole spacing is a common tool for locating data of questionable quality. Thus, pyres includes the option to plot the pseudo-section ( plotutils.plotpseudo) or the raw resistance data (plotutils.plotraw), where plotting of the raw data can incorporate simple quality control tests.
After the ER model has been run, the results are loaded into Python for post-processing analyses using the pyresutils.loadfwdoutput or pyresutils.loadinvoutput functions for forward and inverse model data, respectively. For plotting the ER inversion results, plotutils.plotres or plotutils.triplotres both load and plot the data for ER models with quadrilateral or triangular elements, respectively.
Examples 3.1. Forward–inverse models from R2 documentationPython scripts for all eight surface electrode array examples from the R2 documentation are included with the pyres distribution (Binley ).
These scripts replicate the forward and inverse model input parameters for the R2 distribution examples to demonstrate the capabilities and options in pyres. While all eight example scripts are available in the pyres distribution, only two of the examples will be explained further here.
For more information on the Occams type inversion performed by R2 and ER surveys more generally, the reader is referred to Binley and Binley and Kemna. The first example of using pyres for R2 forward and inverse modeling is to show how pyres implements line topography on a 25 electrode dipole–dipole array (Surface electrode array 2 in R2 documentation; script: R2examples/Surface2dpdp.py). The synthetic earth model for this example contains a 2 m long by 3 m deep 10 Ω m target block 1 m below the surface within 100 Ω m background half space (figure (a)). The land surface elevation increases in the direction of higher line distances with a slope of 0.01 m/m. The inversion converged after two iterations and correctly located the lower resistivity target (figure (b)). However, the inversion overestimated the target resistivity by approximately a factor of three (2.7–3.1) and was quantified using the normalized misfit (figure (c)).
Example of pyres implementing a forward and inverse model using R2 including surface topography. (a) The forward model domain with a 10 Ω m target in a 100 Ω m half space. (b) The inverted ER tomogram from running the forward model. (c) The normalized misfit for the inversion results defined in equation.The second example of using pyres to model ER data with R2 demonstrates how to construct a finite element mesh comprised of triangular elements (Surface electrode array 8 in R2 documentation; script: R2examples/Surface8dpdp.py).
The earth model again contains a 10 Ω m conductor in a 100 Ω m half space, but there is no topography included (figure (a)). The construction of the triangular mesh is achieved by calling meshtools.meshR2 with options detailing the electrode spacing, number of electrodes per line, target information, and discretization options for the domain foreground and background. Topographic information for the land surface can be supplied to the meshing procedure, as in the field data example in section.
Like the previous example, the inversion converged in two iterations and located the conductive target while overestimating its ER. Similar to the.vtk outputs created by R2, pyres can plot the inversion results on the triangular mesh using the plotutils.plottrires function. Example of a forward and inverse model using triangular finite elements with (a) the prescribed half space and target and (b) the inverted tomogram.
Field data inversionField ER measurements can also be inverted and used to construct the input finite element mesh using pyres. As an example, a 56 electrode dipole–dipole survey is loaded and inverted in R2examples/FieldDataExample.py, which was collected near a beach in Rarotonga, Cook Islands (Befus et al). This dataset was chosen as the field example and for later analysis of the sensitivity of the inversion to the error model parameters since it contains three distinct ER targets: (1) a vadose zone, (2) a freshwater lens, and (3) a saltwater wedge. The ability to accurately determine the locations and features of such targets is important for hydrologic interpretation.A mesh with triangular elements is constructed using a characteristic node spacing of one-sixth of the 1.25 m electrode spacing. The land surface topography is loaded using a delimited text file. Unlike the previous synthetic examples, this field survey could contain noisy or erroneous data points. The plotutils.plotraw function can both plot and filter the raw measured resistances, which for this dataset do not contain any poor quality data that exceed 75% misfit in seven point moving kernels.
The inverse model converged after two iterations (figure (a)). The unitless model sensitivity showed the expected decrease in modeled ER values to the measured resistances with depth (figure (b)) and can be used to mask areas of the inverted ER model that contain poorly constrained values (figure (c)). In this example, the inverted data were bilinearly interpolated to a regular grid prior to plotting instead of plotting triangular elements, as was done in the second synthetic model example.
Example of inverted field data after Line L4 in Befus et al. (a) Tomogram of inverted results, (b) model sensitivity data (i.e., diagonal of sensitivity matrix), and (c) transparency applied to inverted tomogram using a minimum sensitivity value of 10 −5 and full sensitivity as 10 −2. Mesh convergenceOne of the goals of ER inversion testing is to maximize the fidelity of the inverted tomogram while using as coarse a mesh as will appropriately model the ER problem.
Too coarse of a mesh could cause the inversion to converge slowly and not solve for resolvable ER features. Whereas, a very fine mesh could resolve all of the resolvable ER features while requiring significant computational resources that were not strictly necessary to resolve those features. For only modeling the forward ER problem, the mesh discretization can also affect the resulting resistance measurements.Given the importance of constructing a mesh that is computationally efficient while still being capable of using the data fully, mesh convergence analyses can be performed to test how well differently constructed meshes perform for inverting a specific dataset. Mesh convergence tests use increasing finer computational meshes to solve the same mathematical problem to establish the mesh resolution required to model the problem ‘sufficiently accurately’. As an example of testing mesh convergence, 105 different inversions of the Surface electrode array 2 forward problem with topography and a buried conductor were performed using increasingly finer meshes, defined by an increasing number of mesh divisions per electrode spacing. For each inversion, the normalized misfit (equation ) between the known earth model and the inverted model were calculated and then integrated over both the whole model domain as well as separately over the elements defining the conductive target (figure ). The integrated misfit over the whole domain quantifies how well both the target and the surrounding medium are resolved using a given mesh, while the misfit for only the target area shows how accurately the inversion predicts the ER of only the target area.
Alone, the conductive target converges to stable integrated misfit values when three mesh divisions per electrode are used, suggesting at least three mesh divisions per electrodes are needed to accurately model the target. However, the whole domain misfit only stabilizes with more than seven mesh divisions per electrode, where eight mesh divisions were surreptitiously used in the original example (figure ). A mesh convergence test for a buried conductor with sloping surface topography (Surface electrode array 2; figure (a)). The integrated misfit decreases as a function of smaller mesh elements for the total solution domain after seven and for the conductive target to about three mesh divisions per electrode. The discretization of the forward model used to create the input data for the inversion also affected the integrated misfit, especially when considering the whole domain.Thus, the mesh discretization is an important control on the inverted ER structure, especially for increasingly coarse meshes. For large ER datasets, this effect of the mesh on the inversion results becomes more important; as more computational resources are required, and the mesh could be constructed too coarsely to accurately model the ER structure in an effort decrease the computational needs. Often, ER inversion software prescribes default mesh spacing based on the electrode configuration, potentially overcoming this potential mesh convergence problem.
The pyres package facilitates conducting mesh convergence tests by operating within a programmatic environment, but mesh convergence tests can and should be performed using other ER inversion software to ensure defensible ER results and interpretations. While not performed in the mesh convergence example, the number of mesh divisions can change in different directions, where the effects of mesh anisotropy could also be investigated.
Alternatively, irregular quadrilateral meshes with non-uniform column and row spaces could be tested using a more constrained mesh discretization (e.g., a R2 ‘type 5’ mesh). Sensitivity test. 2where ER( x, z) i refers to the inverted ER models and m ir are the constant reference models with i = 1, 2 representing two models. While the DOI analysis is useful and convenient to calculate (and is included in pyres: pyresutils.oldenburglidoi), other parameters in the forward and inverse modeling steps can also influence the resulting synthetic data or tomogram, including the electrode configuration, noise in the measurements, and the model parameters (Day-Lewis et al). The pyres package facilitates developing extensive parameter sweeps by looping over parameter ranges and combinations of interest, ultimately running the ER model for each parameter combination. Similarly, for planning ER surveys, the resolvability of resistive or conductive targets can be explored using deterministic or statistical distributions of the target dimensions, line location, and depth along with varying electrode spacing.
3with R the measured resistance values. Thus, a defines the offset error of the data in the units of R Ω, and b is the relative error in R and is unitless. For the sensitivity test, a was log-linearly varied from 10 −4 to 10 Ω, where the maximum measured resistance was 32.0 Ω and the median 0.02 Ω, and b varied from 10 −2 to 0.5 along with a and b equaling zero. Thus, the field data were inverted 132 times to perform the sensitivity test using all combinations of a and b (figures and (a)). When both a and b are zero, R2 uses the standard deviation in R from each measurement and ideally collected with reciprocal measurements and not from repeat measurements. However, only repeat measurements were available from the example field data, so the standard deviation in R was estimated by multiplying the repeat measurement error by the measured R.
It is important to note that not all combinations of a and b will lead to meaningful model convergence (i.e., where the starting model is treated as ‘converged’ due to large a or b) or low root mean square error (rms) values. Each inversion used the same triangular finite element mesh and was allowed up to ten outer iterations to converge.The convergence of the inverse models showed the expected behavior of poor convergence at both high and very low a and b values along with more iterations required for convergence when either a or b was low (figure ). When the error variance model was constrained only by a or b with the other parameter set to zero, the inversion did not converge after ten iterations for a. Convergence behavior for the inverted field data with varying parameters for the error variance model, a and b in equation.Comparing the inverted ER tomograms qualitatively and quantitatively can establish which error variance parameter combination(s) to select for interpreting inversions and aid in deciding which features in the tomogram describe the true subsurface ER structure. Nearly all of the inverted tomograms show a thin resistive layer overlying more conductive material (figure (a)), which has been interpreted as a thin vadose zone overlying saturated sediment (Befus et al). Additionally, towards the beginning or left side of the tomogram, the resistive layer thins and nearly the entire tomogram becomes increasingly conductive for all of the inversions (figure (a)), though the magnitude of these changes depends on the a and b values. This decrease in ER also matches the expected hydrogeologic setting, where the low line distances approach a beach with an expected saltwater–freshwater interface in the shallow subsurface (Befus et al).
These three broadly shared features of the inversions suggest the subsurface can at least be separated into three distinct zones, the boundaries of which depend upon the choice of a and b as well as cut-off ER values. However, many additional features could be interpreted within these broadly defined zones. For example, the wedge-like low ER values on the left side of the tomogram becomes increasingly distinct once a 10 −2.5 Ω would suggest a more diffuse interface. Beneath the shallow resistive layer, the sporadic lobes of higher resistivity represent another example of finer scale features that could be interpreted as the inversion solves for a coarser model with lower a values.
These features may indeed have a physical cause, such as free convection (Stevens et al, Van Dam et al) or evidence of other freshwater–saltwater mixing, but they could also be artifacts of the inversion resulting from overfitting the data. Ultimately, the choice of which ER features are real or not depend upon the choice of the human operator as well as any additional information on the structure of the subsurface, where such ground-truth is crucial for overcoming the non-uniqueness of ER inversions.
Results of the inverted field data showing the (a) inverted tomograms using the same color scheme as figure and (b) the normalized misfit of the inverted models relative to the inversion shown in figure using a = 10 −2 Ω and b = 0.03 (outlined in red). The tomograms include the full foreground domain and are not masked like the tomogram in figure.While the qualitative comparison of ER patterns between inversions can offer some insight into which ER features are real, quantitative comparisons can also indicate the potential variability in the model resolution or range in ER for a known or suspected feature. The choice to save the sensitivity or resolution matrices from the inversion are options in R2, and both provide different information on how well the model can predict the ER throughout the tomogram (Binley ), providing solution-based measures of the inversion performance. Alternatively, the degree to which different inversions diverge from a best-estimate tomogram offers quantitative information in how the error variance parameters change the resulting tomogram (figure (b)).
Using the inversions from the sensitivity test for the field data, the combination of a and b used in the earlier example (figure ) was set as the base model from which all other inversions were compared using equation. These profiles of the normalized misfit highlights the divergent features in the tomograms that arose from changing the error variance model parameters, again showing both the introduction of high ER patches with decreasing a and extensive ER increases with a ≥ 10 −1 Ω (figure (b)).
For this dataset, using 10 −1.5 ≤ a ≤ 10 −2.5 Ω and 10 −2 ≤ b ≤ 0.2 led to relatively small changes in the overall ER pattern and less than ∼50% difference in the ER magnitude from the base tomogram. While it is important to understand how a and b alter the inverted tomograms, many other parameters in R2 will also influence the resulting tomogram and are not independent from the ER structure being surveyed and the quality of the measurements (Day-Lewis et al), where several parameters could be considered in a sensitivity test.
ConclusionsThe pyres package encapsulates the R2 ER modeling software with extensive input and output handling capabilities that expands the utility of R2 as a research tool. Users can readily explore how mesh and model parameters affect resulting ER forward and inverse outputs. These tests can then supply the quantitative foundation for assessing uncertainty in ER models, not restricted to those inverted using R2, and inform their interpretation to enhance the reliability and accuracy of ER surveys in solving scientific and engineering problems. In the future, pyres may also include routines for managing data for the related three-dimensional ER modeling software R3t, and additional utilities for pre- and post-processing ER data may be developed and added to pyres via community contributions.
AcknowledgmentsKMB was supported by the National Science Foundation (EPS-1208909). The pyres package is available from.