NUMERICAL MODELING OF ELECTRICAL RESISTIVITY DATA USING LAGRANGIAN INTERPOLATION

NUMERICAL MODELING OF ELECTRICAL RESISTIVITY DATA USING LAGRANGIAN INTERPOLATION

ABSTRACT This research has examined the numerical modeling of electrical resistivity data using Lagrangian interpolation. The Lagrangian generated a polynomial of minimum degree that represents/mimic the original trend of the data. The reason behind the formulated lagrangian is to reduces the stress of interpretation and improve the smoothness of the field curve. The art of the modeling involves a range [a,b] R when R is the range covering the entire data. The numerical concept is applied to reflect the real trend of the data collected. The Lagrangian interpolation skill was carried out using the MATlab program, and the curves generated showed a good fit that represent the entire trend in the range [a,b]. In respect of this, three VES locations were used as test, to confirm the efficacy of the Lagrangian. All the locations showed that the Lagrangian is efficient with the error values maintained at 2.5%,5.6% and 6.6% compared with the conventional error maintained at 2.5%, 6.4% and 5.8%. The research however concluded that the lagrangian of minimum degree are intelligent to guest the trend of the whole data when considered as a whole at values [a,b] of the generated polynomial.

CHAPTER ONE INTRODUCTION 1.0 BACKGROUND TO THE STUDY Electrical resistivity modeling is an essential tool in the evaluation of trend of electrical resistivity data. It enables to predict the behavior of the electrical resistivity system in response to future stresses due to abstractions. The numerical modeling using lagrangian interpolation, entails the use of a polynomial of minimum order, which represents the trend of the electrical resistivity data. The research data used in testing the efficacy of the lagragian is collected at Alagbaaka GRA, Akure city in Ondo state in Nigeria. Three different locations within the location area is analised and compared with the already interpreted data. The lagragian ensured that closer conceptual model to reality is used, which made the numerical model to be more accurate. The numerical model was simulated based on the conceptual model and utilized code of Matlab7 software as code environments for data input and output management. Once the input data is assessed, conceptual model is developed to represent the field electrical resistivity data system in a simplified version for better understanding of the site conditions. The conceptual model is schematized as a cross section or block diagram which reflects the system behaviors and it leads to determine numerical model characteristics.

(L/2)

(L/2)

Calculated apparent resistivity

Observed apparent resistivity

INVERSION

Final model

Start model

Fig1.0

Block diagram that reflect numerical model characteristics

1.1 STATEMENT PROBLEM The understanding of the electrical resistivity curves is always uncertain. The non-uniqueness of electrical resistivity curves model solution is due to various combinations of data collection, geophysical configuration and geometry factor. Therefore, the introduction of these additional parameters with sufficient number of solute concentration observations in numerical model can

be used to constrain and validate model. Thereby the degree of freedom is reduced in the model, thus improving the reliability of the model using Langragian. 1.2 PURPOSE OF THE STUDY The detailed objection to this study is to model a physical system using lagrangian interpolation to represent the Earth model of the area. Apart from this broad objective, this study also intend to: 

Compare the model curves on the field curves in other to obtain a more fitted Earth curves with minimum error, and high reliability.



Design a system that can be used to interpret other basement complex environment.

1.3 THE SIGNIFICANCE OF THE STUDY The research will be useful in analyzing different types of curves in an environment that is subjected to conditions stated in the Earth model.

1.4 DEFINITION OF TERMS

i.

Array: Arrangement of configuration

ii.

Curves: is a graphical representation of field data

iii.

Conventional curve: is a curve matching system

iv.

Data: is a set of values of qualitative or quantitative variables

v.

ERT: Electrical resistivity topology

vi.

Model curve: is a curve that trend the conventional curve

vii.

Nodes: a point in a network where lines cross or branch

viii.

VES: Vertical Electrical Sounding

CHAPTER TWO LITERATURE REVIEW AND THEORETICAL FRAMEWORK INTRODUCTION 2.0 ELECTRICAL RESISTIVITY METHOD The electrical resistivity method uses a low frequency power source to create an electric field by injecting current into the subsurface via conductive electrodes. The voltage potential is measured at other electrode locations on the ground surface. Modern, multi‐channel acquisition systems allow for various voltage measurements to be made for one current injection and current can be injected using any combination of the available electrodes. Common measurements require four electrodes, two used as current source and sink (A and B), and two used to measure the potential voltage drop (M and N, Figure 2.1). Solid lines in the subsurface display the current flow path as well as dashed lines displaying the equipotential

Fig 2.0: Illustration of the Electrical Method showing current being injected from A to B while voltage potential is measured at electrodes M and N. Based on a combination of measurements, the electrical resistivity method senses the distribution of subsurface electrical resistivity, or its reciprocal, electrical conductivity. Resistivity is an intrinsic material property defined as how strongly a given material opposes the flow of electric current. It is a property that can cover several orders of magnitude. The electrical resistivity method has been used in exploration of the subsurface for many years. The simplicity of the equipment, the low cost of carrying out the survey and the abundance of interpretation methods make it a popular alternative to drilling and testing. Although there are many variations of resistivity technique, the basic procedure is to establish a subsurface current distribution by injecting current into the ground between two electrodes.

A series of potential differences are measured between pairs of potential electrodes in a line or grid, and are then interpreted on yield information about the electrical conductivity beneath the study area.

2.1 APPLICATION The resistivity techniques has been extensively in the detection and mapping of groundwater contamination (American society of civil Engineers,1961). Considerable work has also been done by apply the technique to the evaluation of groundwater resources and in the quantitative evaluation of aquifer parameters that control groundwater flow. Resistivity measurements have also been used in the geostatistical extrapolation of aquifer parameter away from well control. The use of geoelectric methods in the monitoring of enhanced oil recovery (EOR). The use of resistivity survey in mining exploration has been less widespread. This is primarily due to small influence that disseminated mineralization often has on resistivity, as well as strong competition from the induced polarization method. The techniques are also see some use in the exploration for hydrocarbons. 2.2 FIELD PROCEDURE Numerous ways of arranging the current and potential electrodes in a resistivity survey are possible. Both the geometric configuration of the electrodes (i.e. the array used) and the relative spacing between electrodes can vary from survey to survey. Because both the depth penetration and vertical resolution depend on the particular arrangement of electrodes used, a proper choice

of array is important to the success of the experiment. This choice depend on the target (such as depth of burial and resistivity contrast), as well as on the complexity of the surrounding medium. The most common way of acquiring resistivity data is to make use of one or more traverse lines. To image lateral changes in the subsurface, the center of the array is moved along one of the traverse lines while keeping the electrode separations fixed. Arrays that are frequently encountered include the Schlumberger, pole-dipole, dipole-dipole array. Expanding the electrode array about a fixed point to image successively deeper regions of the Earth is referred to as depth surrounding. Combining two techniques into a depth sounding/profiling survey is also possible. 2.3 RELATIONSHIP BETWEEN ELECTRICAL RESISTIVITY AND PARAMETER The close relationship between electrical current flow and porous media transport has stimulated considerable interest in the application of DC resistivity methods of groundwater resource and contamination problems. Much attention in recent years as focused on the relationships between electrical resistivity and parameter governing the porous media flow. The starting point for most studies is to relate the bulk electrical resistivity b to the resistivity of the pore fluid f using a constant formation F, where F= b/ f …………………………………………2.1 Assuming the formation factor is indeed constant, this allows one to determine the ionic content of the pore fluid, and hence the water quality, once the bulk resistivity has been estimated. The assumption of a constant formation factor, as pointed out by ( Spitzer,K. 1998), will be invalid for problems where matrix conduction is significant. In these cases, more complicated formulae

must be considered. For problems where the porous medium is non-homogeneous, theoretical relations like Archie’s law F=av-m ……………………………………………………. …………….2.2 Can be used to relate the formation factor porosity v. combining equation 1 and 2, obtain V=logm(a f/ b) ……………………………..2..3 Equation 3,can be used to determine the porosity of the medium, given an estimate of its bulk electrical resistivity.

2.4 INTERPRETATION The first step in a typical interpretation of resistivity data is to convert the measure potentials into apparent resistivity – the apparent resistivity,

a,

being the resistivity of a homogenous earth

corresponding to the particular potential measurement

. For the pole – pole experiment

the apparent resistivity is given by

a

=2π/I. √

.

………………………….2.4

Apparent resistivity data obtained from profiling surveys are usually displayed as pseudosections or contoured plan maps and interpreted quantitatively. Analytic solutions for spheres, vertical dikes and other simple earth structures are often used to aid in the interpretation. To model more complicated structure, 2D and 3D forward modeling techniques are often used. Recent work on the forward modeling DC resistivity has concentrated on the use of standard

numerical procedures, including the Langrange interpolation, finite difference method (Zhao,S. and Yedlin, M.J. 1996) and boundary element. A number of attempts have been made to solve the inverse problem of 2D and 3D conductivity distributions. Examined the inversion of resistivity and IP data over 2D structures using a data base of synthetic data generated for different models. Because of the need to restrict the size of the data base, only 9 parameters were actually used to describe the model. These include the thickness of the overburden layer, width and height of a buried prism, and resistivity and polarizability of the overburden layer, prism and host. The 3D inversion scheme based on the use of alpha centers. Zhang,j., T.R 1995 used a least-squares approach to solve the 2D inverse problem for resistivity and IP data set. Although the finite difference algorithm they used in the inversion permits a highly detail model, they chose to solve the conductivity of only a small number of sub-regions. The algorithm for inverting the schlumberger data over a 2D earth. Again, the conductivity of only a small number of sub-regions was solved for the inversion. Ellis, R.G Oldenburg,D.W. (1994), using a similar strategy, examined the inversion of cross borehole data for 2D conductivity structures. Because of the computation problems encountered solving the 3D problem; the conductivities of only a small number of cells were solved for in the inversion. Although, the use of inversion techniques for the interpretation of 2D and 3D problem is highly desirable, good results have limited by small number of parameters that can be solved for. 2.5 ERT IMAGING FORWARD MODELING AND INVERSION Forward modeling is calculating ERT data based on assumptions of resistivity structure (spatial) and measurement locations. The inputs for forward modeling are defined in the model space. The

position and magnitude of the resistivity structure, the size of the grid cells, the approximation method, and the locations of the sources and receivers are factors that influence the forward modeling approximation of ERT data. Forward modeling can be used as a pre‐investigation tool to determine whether a geologic target can be resolved using the ERT method. Forward modeling aids in ERT survey design such as electrode spacing, measurement geometry, and optimized acquisition sequences.

Inversion is the mathematical process that allows us to estimate a true resistivity distribution based on measured data. This is in contrast to forward modeling that allows us to estimate data assuming we know the subsurface resistivity distribution. Inverse problems are non‐unique because there is a desire to find a larger number of model parameters from a limited number of observed data points. Inversion is an iterative mathematical process whereby the subsurface structure is adjusted until the forward modeled data (dcalc) matches the field observed data (dobs). The resulting inverse model is a non‐unique solution that honors the observed data to within some measure of fit.

2.6 DATA SPACE In practice we measure the voltage drop in a circuit as well as an injected current, which results in an estimate of the resistance of the circuit. ERT data is referred to as apparent resistivity, app, because it is a measured resistance 𝑅 =

in Ohms scaled by a geometric factor K (m)

which accounts for the measurement geometry. The geometric factor (K) is based on the assumption that the earth is a homogeneous, isotropic half‐space and that the resulting

equipotential surfaces around current source and sink form symmetric hemispheres (Broyden,C.G 1972). (R app = 𝐾 The geometric factor does not account for topographic variability or heterogeneity in the subsurface but nevertheless “apparent” resistivity in a common way to scale measured resistances.

2.7 DC FORWARD MODELING The governing equation that relates the electric potential field from a unit input current dipole is provided.

‐∇ ·σ(x,y,z)∇ (x,y,z)=I(δ(r‐r+)‐δ(r‐r‐))…………………………………………………2.6

Variables from equation above are defined as: ∇ · ‐ divergence σ ‐ the conductivity structure of the medium ∇ ‐ gradient Φ ‐ is the electric potential field I ‐ is the input current from a dipole r+ and r− are the locations of the positive and negative current sources δ (r‐r+) is the Dirac delta function, centered at the current source location

The equations are implemented in Matlab as described in FW2_5D: A MATLAB 2.5‐D electrical resistivity forward modeling code (Zhang,J.,Mackie,R.L., and Madden, T.R. 1995). The program goes through a series of operations as well as transforming to and from the Fourier domain to estimate the electric potential field given a resistivity structure and the locations of current sources and sinks. For 2.5D, the solution is three‐dimensional but the third dimension is invariant.

(DS( )G)u=A( )u = q………………….2.7

In matrix form, D and G are divergence and gradient matrix operators, S( ) is the matrix containing the conductivity structure, u is a vector containing the potentials, A( ) is the combined forward operator, and q is a vector containing the location of the current sources and sinks.

dest = Qu = QA1( )q…………………..2.8 The forward modeled data (dest) is calculated by multiplying the potential field vector (u) by a projection matrix Q, which selects a subset of the potential field based on the positions where the voltage potential is measured in the survey. For a 2D line profile we can only make voltage measurements at the surface.

Boundary Conditions The forward modeling algorithm explicitly models the ground surface/air boundary by assigning bulk resistivity to the cells in the ground as well as air resistivity value (1x106 Ωm) for the cells above the ground. The edges of the model space are no‐flux boundaries although the model space is padded with extra boundary cells.

Grid Setup and Dimensions The boundaries of the model space extend from ‐370 m to 370 m centered around zero while the top of the model is zero and the bottom edge is 55 m using a positive downward convention. The 72‐electrode array extends from ‐141 m to 141 m laterally and the deepest point of the line is located in the base of the pond at 6 m. To reduce numerical error in the forward model, the grid is finely discretized near the electrode locations. The size of the grid, including the pad is 1,294 cells wide by 197 cells deep. The cells are variable in size but the smallest dimensions, located near the electrodes are 0.25 m wide by 0.05 m deep.

2.8 INVERSION ERT inversion is the process of reconstructing a resistivity structure based on acquired data. Solutions to inverse problems are non‐unique for multiple reasons including ill‐posedness, data noise, and in particular the larger number of model parameters, m (100,000’s), from a limited number of observation points, dobs (100’s ‐1000’s). All ERT inverse algorithms are based on a measure of fit or comparison between the observed data (dobs) and calculated data (dcalc) as well as perceived structure (mref) versus the inverted structure (m). The perceived structure is

based on “a priori” information you may have about the subsurface conditions. Mathematically, where Wd is a dataweighting matrix (Zhang,J.,Mackie,R.L., and Madden, T.R. 1995).

Data Norm = ⎢ ⎢ Wd (dcalc‐dobs) ……………………………2.9

The model misfit is given by equation 3.2 where Wm is a model‐weighting matrix (Zhang,J.,Mackie,R.L., and Madden, T.R. 1995).

Model Norm = ⎢ ⎢ Wm (m‐mref) ⎢ ⎢ 2 …………………………………………..2.10

The objective function is a combination of data fitting and “a priori” model fitting. The objective function uses a combination of zero and first order Tikhonov regularization (Zhang,J.,Mackie,R.L., and Madden, T.R. 1995). Solving the inverse problem is an iterative process whereby the resistivity structure is updated until the objective function,

, is

minimized (Zhang,J.,Mackie,R.L., and Madden, T.R. 1995).

⎢ ⎢ Wd (dcalc‐dobs) 2 +

Wm (m‐

The trade‐off parameter, is a user‐defined value, which balances the relative importance of data misfit vs. realistic model structure. In this research we set the trade‐off parameter to a level consistent with previous inversions completed on the HSRP ERT dataset. It may require several iterations of the resistivity structure before the calculated data fits well enough to the observed data. The data and model norms are the respective components of the objective function. The

model norm is measure that represents the difference between the inverted resistivity structure (m) and the reference structure (mref). The forward and inversion algorithms are solved using a finite volume method on a regular grid in Matlab. The inversion program is a 2.5D version of the open source code (Zhang,J.,Mackie,R.L., and Madden, T.R. 1995).

The data‐weighting matrix is commonly chosen to be a diagonal matrix of 1/𝜎, where 𝜎 is the standard deviation of the measurements (dobs). This term represents the error model for our data. Repeat measurements can be used to quantify individual measurement errors.

Wd = 1/(3%* dobs) …………………………………..2.12

In the absence of repeat measurements, a common assumption is that measurements follow a Gaussian distribution and the standard deviation of a single measurement is 3% of its absolute value (Archie, G.E,1942). With this assumption the data‐weighting matrix (Wd) normalizes the measurements by their respective standard deviation so that measurements of varying magnitudes have a similar influence on the inversion result. Without weights the least‐squares (L2) norm would be biased to fitting the larger magnitude measurements.

3

The model‐weighting matrix (Wm) is comprised of three terms that govern the smoothness and stability the solution (Zhang,J.,Mackie,R.L., and Madden, T.R. 1995). Wm is often referred to as a roughness or smoothness parameter as it can effect how sharply or smoothly

‐4) and an identity matrix (I) for numeric stability. Wx and Wz represent the gradients (first‐order Tikhonov) in the x and z directions, respectively (Zhang,J.,Mackie,R.L., and Madden, T.R. 1995). As the trade‐off parameter , the objective function is minimizing the data misfit regardless of how physically unrealistic or blocky the resulting resistivity structure is (i.e. the first term in Equation 3.3 becomes dominant). If it is very large, then the resistivity structure will be smoother but may have a large misfit from the observed data (i.e. the second term in Equation 3.3 becomes dominant). Equations provided in this section are modified from RESINVM3D.

2.9 SIGNIFICANCE OF MODEL INPUT PARAMETERS We have shown in Chapter 2 that the pond stage dominates the ERT signal when compared to the same synthetic resistivity structure with no stage. In this section we examine the effect that a slight variation in pond stage can have on ERT data. The question is, “What effect do errors in the position and/or timing of pond stage data and surveyed topography have on the approximation of the resistivity structure?” It is well known that electrode positional errors will propagate through inaccurate forward modeling of the data (Oldenborger et al, 2005). The purpose of this exercise is to quantify the forward modeling error in this scenario where a number of our electrode locations are in close proximity to the pond stage interface.

We demonstrate the effect that stage and topography errors can have give three scenarios with homogeneous, dry sand subsurface and a ponded stage layer. Model scenario A is the baseline scenario (i.e.: no error in either topography or stage), while model scenarios B and C simulate practical measurement error in pond stage and topography, respectively.

Table 2.1: Modeling Error Scenarios Model Input

Description

Data Output

Resistivity Structure A

Baseline synthetic model with a dry sand sub and a ponded water layer

dcalcA

Resistivity Structure B

The estimated pond stage is moved 20 Cm higher than it actually is in A

dcalcB

Resistivity Structure C

A single electrode and its topography is moved 20 cm lower than it actually is in A

dcalc

2.10 POSITIONAL ERROR For the 2D ERT surface line we can expect positional errors on the order of +/‐ 10 to 20cm in the surveyed topography or perhaps larger. Positional errors are not uncommon when using a total station for surveying and there can be a bias in the error if the total station has not been leveled accurately. The bias is such that the positional error increases as your survey target is further away from the total station. Figure 2.1. illustrates how a biased error in electrode positions may propagate compared to their true locations with the total station positioned on the left flank of the pond.

Fig 2.1:

Diagram showing a potential error scenario in electrode survey positions.

Random positional errors due to operator inexperience with the total station survey technique may also be an issue. A comparison between the true topography and random estimate of topography is shown in Figure 2.2.

Figure 2.2: Diagram displaying a potential random error scenario in surveyed electrode positions.

2.11 MODELLING This section presents a technique for the efficient numerical computation of the electrical potential with the finite element method in three dimensions and arbitrary topography. The crucial innovation is, firstly, the incorporation of unstructured tetrahedral meshes which allow efficient local mesh refinement and most flexible description of arbitrary model geometry. Secondly, considerably more accurate results are achieved by the implementation of quadratic shape functions. Exploiting a secondary potential approach, meshes are downsized significantly in comparison with highly refined meshes for total potential calculation. However, the latter is necessary for the determination of the required primary potential in arbitrary model domains. To start with, homogeneous models with different geometries at the surface and subsurface are simulated to quantify their influence. This results in a so-called geometry effect, which is not only a side effect but may be responsible for serious misinterpretations. Moreover, it represents the basis for treating heterogeneous conductivity models with the secondary potential approach, which is especially promising for the inverse problem. It is addressed how the resulting system of equations is solved most efficiently using modern multi-frontal direct solvers in combination

with reordering strategies or rather traditional preconditioned conjugate gradient methods depending on the size of the problem. Furthermore, a reciprocity approach to estimate modelling errors is presented, and it is investigated to which degree the model discretization has to be refined to yield sufficiently accurate results.calculation of the electric field started in the late 1960s using the techniques of integral equations and finite differences (Bing, Z. and Greenhalgh, S.A. 2001). A special variant based on integral equations is the boundary element method (Farqharson, C.G and Oldenburg,1998). Finite difference (FD) calculations in three dimensions. They were the methods of choice throughout the 1980s and 1990s and were refined several times. Zhang et al. (1995), e.g., presented improved boundary conditions for a more accurate potential approximation and Spitzer (1995), e.g., introduced efficient preconditioned conjugate gradient solvers to decrease execution time. The finite difference method is restricted to orthogonal grids which limit its ability to reproduce non-orthogonal geometries. In recent years, more and more finite element (FE) approaches appeared which are generally not subject to these drawbacks. FE formulations of the direct current (DC) resistivity problem are described in Li and Spitzer (2002); Zhou and Greenhalgh (2001). Simulations with surface topography were presented in two dimensions (2D) (Zhang,J.,Mackie,R.L., and Madden, T.R. 1995) and three dimensions (3D), Yi et al. (2001). However, the presented algorithms mainly work with block-oriented (structured) discretizations using hexahedral or tetrahedral elements and, thus, do not exploit the full power of the FE approach. Sasaki (1994) and Zhou and Greenhalgh (2001) use tetrahedral grids but since they are derived from bricks the approach is still blockoriented. The grids in (a) and (b) are of type regular (structured). The orthogonal hexahedral grid (a) is the one that is furthermost restricted with respect to geometrical flexibility and local refinement. The local refinement of (b), a non

orthogonal regular hexahedral type, is still awkward and inefficient, but its geometrical adjustability is already increased. In this work the C++ class library GIMLI has been developed, which provide the solution of finite element problems and applies the non-commercial mesh generator TetGen Si (2004) for the generation of unstructured tetrahedral meshes. Therefore two subjects can be addressed. Firstly, any model topography and electrode layout can be described flexibly, e.g., arbitrary electrode positions can be represented by prefixed nodes within the mesh. Secondly, the node distribution can be controlled on demand by refining the mesh locally in the vicinity of electrodes or at strong conductivity contrasts, i.e., where strong gradients of the simulated potential require enhanced accuracy. In return, coarse grids toward the boundaries are sufficient to approximate the rather smooth fields. As a result of any discretization technique, a sparse system of equations has to be solved for each current location. Spitzer and Wurmstich (1999) gave an overview on non-stationary iterative equation solvers for DC resistivity problems. They suggest conjugate gradient techniques with preconditioners. One very efficient, albeit memory consuming technique is obtained by an incomplete Cholesky decomposition as applied to DC modelling by Zhang,J.,Mackie,R.L., and Madden, T.R. 1995 and Li and Spitzer (2002). In contrast to these iterative methods direct equation solvers have made substantial progress in recent years. Multi-frontal decomposition methods going back to suit best for problems with many right-hand sides, since the decomposition is done only once. The use of reordering techniques helps to limit the memory requirements for direct methods as well as for the incomplete Cholesky preconditioner. Accuracy in DC resistivity modelling is enhanced not only by grid refinement strategies but also by employing higher order base functions in the FE approach. It can be shown that a

combination of both, mesh refinement and quadratic shape functions, yields most efficient solutions. First an introduction to finite element modelling, unstructured mesh generation and refinement is given. The homogeneous half-space is used to investigate how refinement strategies improve the results for the calculation of the total potential. As a continuation, a conducting half-sphere is used to demonstrate the geometric flexibility of the presented approach and to apply the secondary potential modelling.With the help of two examples it will be shown how complicated subsurface and surface geometry can be involved. A geometry effect is defined to appraise the influence of any topography on DC resistivity measurements. Finally, the computational aspects for the solution of the system of equations are discussed.

2.12 CALCULATION OF THE TOTAL POTENTIAL The boundary value problem for the direct current (DC) resistivity forward problem is given by the equation of continuity with mixed boundary conditions:

∇ 𝜎∇

∇

………………………………………2.14

………………………………… 2.15

where 𝜎(x;y; z) represents a given conductivity distribution in the ground, j the source current density and u is the sought electrical potential. The electrical resistivity , which is the most

widely used material parameter in DC geophysics, equals the reciprocal of 𝜎. The classical geoelectrical source term for the current density j reads

……………………………………… 2.16

∇

and assumes a current I being injected through a point-like electrode situated at rs = (xs,ys, zs). The Dirac’s delta d forces a singular current density at the source position. The formulation (2.15) allows point sources at arbitrary positions inside the domain as well as sources on the boundary of the modelling domain. However, the expression does not take into account the physical presence of finite electrodes that may influence practical measurements. The type of boundary condition for the geoelectrical problem in equation (2.16) is specialized by . At the earth’s surface δ

(i.e. air=material interface) homogeneous Neumann conditions

( = 0) are applied to avoid current flow through the boundary along the outward normal n.

To allow an asymptotic decrease of the potential on an artificial subsurface boundary δ

,

impose a = (n. r)/r2 for a current source at the earth’s surface, whereas r =|r| is the radial distance between the current source and the corresponding boundary. A more generally form for a buried current source inside the domain is given by (Zhou and Greenhalgh,2001)

= _______________ …………………………………… 2.17 (

(r+r’))

with

=|

where

denotes the image position of r regarding the earth’s surface. This type of

boundary condition is sometimes referred to as Robin type (Herwanger et al., 2004). Note that on pure Neumann domains, e.g.

= 0 on the whole boundary as in modelling tanks, there are two

additional requirements ensuring a unique solution. First, the conservation of Charge (total current flow on the whole boundary δΩ is zero), which can simply be achieved by using a reference electrode or dipole current pattern. Secondly, since only derivatives of u are present on the boundary, the solution of equation (3.1) becomes ambiguous and hence the system of equations becomes singular. This can be avoided by choosing an arbitrary reference point on the boundary with the additional condition u = 0. The Dirac delta function d in the source term of the DC equation leads to infinite potential gradients at the source position rs. Typically, this singularity is responsible for very poor numerical approximations, particularly close to the electrode positions. The total potential (TP) u is split up into a primary and a secondary part, u = up +us. The primary potential up is a known electrical field for a given background conductivity structure sp and satisfies the source term: ∇ 𝜎 ∇

= -∇

Usually the analytical

solution for the homogeneous half-space is used but also more complex background structures (vertical dike or layered models) can be applied (Li and Spitzer, 2002).

The boundary value problem for the secondary potential (SP) reads:

∇ 𝜎∇

∇( 𝜎

𝜎 ∇

)

………………2.18a

+

=

+

I, ………………2.18b

where the singular current density j is vanished. However, the integral of equation contains the gradient of the primary potential up. The precondition for it being regular is that the conductivity in the direct vicinity of the electrode equals 𝜎p, as already mentioned by Zhao and Yedlin (1996). To determine sp, the local conductivity at the electrodes has to be chosen and not the mean conductivity . Secondary sources appear where the conductivity deviates from 𝜎p.

Both left-hand side and right-hand side differential operators are identical to that of the total potential. By approximating the operators by matrices the system of equations can be written as:

us =

up. ………………………………2.19

In order to avoid the assembling of

for each source conductivity , the linearity of A is

exploited by

up = (

Thus only one matrix to

-

)up =

up𝜎 -

up. …………………… 2.20

for a homogeneous conductivity of 1 has to be created additionally

for the whole forward calculation with many sources and heterogeneous conductivity

distribution.

2.13 MESH GENERATION AND REFINEMENT An initial partitioning of the modelling domain often incorporates the given electrode layout as fixed nodes. However, an additional refinement has to be applied to obtain accurate forward calculations due to the singular potential at the electrodes. The existing approaches use block oriented grids for both FD and FE calculations. To minimize the error of the forward calculation, grid lines between the electrodes are introduced. Using between two and four additional grid lines all calculations end up in around 4% relative error for a pole-pole configuration. For a successive superposition of configurations with large geometric factors (such schlumberger) a further refinement is required since the relative error is amplified by the geometric factor. However, since refining block-oriented grids always works in a global way, the number of nodes increases rapidly. One main advantage of unstructured meshes is the facilitation of refining grids within distinct regions. Thus in regions of varying potential gradients (close to electrodes) the mesh can be chosen very fine whereas the cell sizes grow toward the boundaries of the modelling domain. It can be distinguished between a-posteriori and a-priori refinement. For the former, the discretization depends on an error estimation procedure in the solution process, whereas for the latter, the information is introduced in advance. Since the critical regions are known, the a priori type is chosen to enforce a locally fine mesh by introducing additional supporting nodes to the mesh generation process.

Quality measures In addition to the local node density, the approximation quality depends on the cell size growth (or prolongation) factor. This can mainly be controlled by the ratio of the tetrahedral edge lengths and the radius of the circumscribing sphere. Since sliver-shaped elements yield poor approximation properties, so-called mesh generators try to minimize the radius-edge ratio. The maximum ratio throughout the mesh can be used as a global mesh quality control.

The applied mesh-generator tries to force all radius-to-edge ratios below a certain quality constraint. For all subsequent meshes a radius-to-edge ratio of 1.2 is chosen, which, due to experience, provides sufficiently accurate results. To create an unstructured mesh the domain has to be defined by points, polygons or faces, the so-called piecewise linear complex (PLC). It includes the geometry of the domain and the electrode locations as well as the boundary. By introducing nodes and creating tetrahedral elements, the mesh generator creates the mesh of the desired quality.In the following study, local and global refinement techniques are investigated to find a best trade-off between accuracy and computational effort. In addition to the spatial mesh refinement (h-method), the use of higher order shape functions (p-method) was taken into account.

2.14 DISCRETATION FOR CALCULATING THE TOTAL POTENTIAL A uniform earth with a flat surface boundary offers the most simple analytical solution for a homogeneous resistivity distribution

= 1/𝜎 = 1

with:

U = I /2

………………………2.21

at radial distance r(x;y; z) from source location rs(xs;ys; zs). A straight line is defined containing 21 electrode nodes with 1m spacing. The model boundaries are placed at 5 km around the origin in order to minimize the effects of the boundary conditions. The first mesh created by the mesh generator TetGen obtains 2047 nodes.

Discretization for calculating the secondary potential In order to show the flexibility of unstructured meshes a spherical anomaly was chosen, which can hardly be discretized by block-oriented grids. The analytical solution for a conducting sphere is known ((Li and Spitzer, 2002)). However, instead of a sphere in full-space, a half-sphere at the upper boundary of a half-space is considered, on which electrodes are placed along a profile line. Since the half-space boundary is an axis of symmetry, the solution can easily be obtained by doubling the calculated potential values. The secondary potential is expected to be rather smooth so that the first calculation may start with a quite coarse mesh.A halfsphere with a radius of 2.25m is placed at the origin. A line of 21 electrodes from x = �5m to x = 5m with a spacing of 0.5m is defined using fixed nodes. Note that the electrodes can be placed independently of the nodes in general, however, defining fixed nodes in this example avoids possible interpolation issues while collecting the resulting potentials. The model boundary is also chosen away from the origin in each direction (1000 m). Using a radius-to-edge ratio of 1.2 a mesh of 1769 nodes has been generated on which the secondary potential is simulated. Defining the resistivity of half-sphere to 1Wm, whereas the half-space has

= 10Ωm. The source electrode is set at rs = (-4; 0;0) and the primary potential is calculated with equation (2.21) assuming a background resistivity of 10Ωm.

Consequently, sources for the secondary potential occur only within the half-sphere. The calculated potentials at the remaining electrodes are transformed into apparent resistivities using pole-pole geometric factors.

2.15 MODELLING GEOMETRY EFFECTS In general, the measured electrical impedance

is transformed into the apparent resistivity

by means of the geometric factor k. The latter is chosen such that the apparent resistivity equals the true resistivity in case of a homogeneous distribution. The geometric factor clearly depends on both, the electrode layout and the surface geometry. If the topography is non-trivial, k is unknown and can only be assessed numerically. Usually an approximation ka is based on the half-space potentials for surface measurements. For buried current sources the analytical expression can easily be adapted by introducing an image source position r0s using the earth’s surface as a mirror:

=

=

( +

). ……………………….2.22

It is important for any measurement, to distinguish between effects of the subsurface conductivity distribution and artefacts caused by using wrong geometric factors. Therefore, a geometry effect t is defined as the ratio of the voltage differences for the given geometry u and

the approximation ua. From the equality of the apparent resistivities the expression for the geometry effect yields:

t=

=

…………………………2.23

If any surface topography is present, a homogeneous resistivity of

= 1Ωm is assigned to the

model and the potential u is calculated numerically. Since a has to equal

the numerical

geometric factor can be obtained by k = I / u. To appraise geometry effects, t can be plotted for each datum. With a value of t = 1 the measurement is not affected by topography. Values of t > 1 refer to increased apparent resistivities whereas values of t < 1 indicate a decrease.

2.16 COMPUTATIONAL ASPECTS The meshes of the previous sections are now being reviewed with respect to the effort of solving the system of equations. For iterative methods, a stopping criterion has to be defined depending on the accuracy of the solution. Experience shows a relative residual value of

is sufficient

for practical purposes and was used in all computations. Proceeding in the order of the mesh size and starting with the smallest example, the conducting sphere with 1 769 nodes and 21 right-hand side vectors. The stiffness matrix contains 24 427 non-zero entries which is about 14 per row. Due to the small size, a direct equation solver is appropriate. Once the factorization is done the back-substitutions for the individual sources are obtained almost instantly. The Cholesky factor is computed in 1.85 s and contains 750 238 nonzero elements which need about 11.5MB memory if 12 bytes/entry are considered. All methods

of reordering can reduce the nnz drastically, also the computing time is decreased by several decades. AMD performs best with a reduction factor of about 6. The locally refined homogeneous half-space with quadratic shape functions (38 533 nodes) still represents a moderate mesh with a stiffness matrix A of 542 459 non-zero elements. Without reordering, the Cholesky factor obtains 6.84

elements corresponding to 780MB RAM

which almost exceeds the memory limits of a current standard PC (1 GB). However, with AMD reordering the allocation can be reduced to nnz = 1.07

or 120 MB. The AMD reordering

technique was used for all examples because of its great benefit.

The mining gallery (128 169 nodes, 1.8 Mio nnz’s in A, 50 electrodes) represents the next larger problem. The (reordered) Cholesky factor obtains 6:38.

nnz’s or 730 MB. On the present PC

with 1GB memory, 200 000 nodes represent the upper limit for direct equation solvers. The algorithm implemented in the CHOLMOD package needs about 3 s for the symbolic analysis and 52 s for the multi-frontal factorization. The following back substitutions are carried out in 111 s which is about 2 s for each electrode. Hence, conjugate gradient (library TAUCS) is chosen for solving the system of equations. Nevertheless, it is worth to carefully choose the preconditioner since the system is solved for each of the 16 electrodes. The incomplete Cholesky preconditioning, as used by Li and Spitzer (2002), was selected and both variants were examined with different thresholds. As for the complete factorization, the nnz can be drastically reduced by AMD reordering. To summarize, the numerical computations must always find a trade-off between the two resources calculation speed and memory. Mesh size and available computer memory determine the runtime of the solver. Direct methods should be used whenever it possible. Conjugate

gradient methods with incomplete Cholesky preconditioners are the alternative in most cases. The threshold value of the incomplete factorization must be chosen carefully to meet memory limits and to minimize calculation time. For very large meshes, the SSOR preconditioner is of advantage since it allows preconditioning without additional storage memory (Spitzer, 1995). However, convergence becomes slower.

CHAPTER THREE RESEARCH METHOLOGY 3.0 INTRODUCTION The previous chapter discusses the review of different researchers on numerical modeling, this current chapter intends to look at model conception and methodology adopted in arriving at the signatures of the field curves and tries to compare it with curves generated by lagragian interpolation, this will enable us to know if it is a good representation of the of the field curve

3.1 MODEL CONCEPTION The aim is to generate a polynomial of minimum order that represents a curve that approximates other curves and reduce the error to the nearest minimum. To do this we need to: i.

We shall use the linear and quadratic interpolation formula

=

…………………………….…………………………….. 3.1

Re-arrange: y=

(

)+

If we consider two points. (

(

) ……………………………………………....3.2

) and (

) in the plane.

By using the langrange and interpolation basis and interpolation formula, In the same way these is a unique parabola through three points

,

in the plane. One way of

finding the polynomial is to use linear interpolation between two straight line functions similar to equation3. 2 above. This yields curves ponding symmetric form:

Y=

+

+

…………………………...3.3

Generally these two polynomial can be written as:

∏

Where ( ) = These basic polynomial have property that ( ) =

Therefore the langrange interpolation polynomial given by

∑

For a system of linear equations set up as vander munde matrix:

P(n) = =

+ …………………………

+

+ ………………………+

+

[

][

…………………………3.6 +

………………………..3.7

] = [ ] ……………………………………………………..…3.8

3.2 ERROR IN LAGRANGE Assuming function f is sufficiently smooth, we may deduce, by repeated application of Rolle’s theorem to an appropriate function, that the error in the Lagrange interpolation formula is given by

𝑅 (f:x) = f(x) -

…………………………………………………..3.9

(x) =

Where p is the Lagrange interpolation polynomial

∏

And x and the mean value point

lie in the interval [

remains valid for extrapolation with

] spanned by the nodes. Equation 3.10

in the interval by x and the nodes.

3.3 ERROR BAND FOR POLYNOMIAL INTERPOLATION Since

in equation 3.15 depends on x, it follows that we cannot usually find the error exactly, if

we could, then we could also evaluate f(x) exactly. The merit of this langrange remainder is

therefore that we can use it to obtain error bounds, if we can find a bound on the appropriate derivative of f on [a,b]. Thus, if

|

|

M for every x

[a,b], it follows that

| ……………………………………………………………3.11

|

In order to obtain the best uniform bound for this error over the whole interval, it is therefore desirable to choice the nodes, so that the maximum (absolute) value of

Formula:

=

–

(x)

………………………………………………3.12

Equation 3.15 is clearly symmetric in the nodes. This property carries over to divided differences of arbitrary order comparison of the langrangse and divided difference formulas yield:

[

]∑

∑

Applying Newton’s Divided difference formula: The Linear interpolation polynomial agreeing with f at [ ]

[

] …………………………………………………..3.14

The quadratic interpolation f at three point [ ]

[

is given by:

]

is [

] ………..…… 3.15

These are special cases of Newton’s divided difference interpolation formula:

[ ]

[

[

]

[

]

]

The error for this formula is given by: [ [

]

]

3.4 DIVIDED DIFFERENCES AND DERIVATIVES By considering the error term in the langrange and Newton divided difference formula, we can obtain the following relation between divided difference and derivatives of f:

[

]

, where

is some point in the interval spanned by

. for the special case k=1, this is just a restatement of the mean value theorem. Using the limit in equation 3.26 as

leads to the interpolation of divided

differences for repeated nodes: [

]

…………………………………………………...3.18

over [a,b] is minimized. For the interval [-1,1], this is achieved by choosing the chebyshev nodes which are at (

), (k=0,1,………..,N) ………………………………3.19

The optimal nodes for a general interval [a,b] are at ( ) [

]

3.5 DIVIDED DIFFERENCE For any of the nodes

, the zeroth-order difference is defined by f[ ]=f

difference at the nodes

f[

]=

[ ]

[

.The first divided

is then defined to be:

]

(

=

)

…………..……………………………3.20

and the recursive nature of the definition extends to the general k-th divided difference at nodes .

f[

]=

[

]

[

]

=

………………….……………………3.21

Note that the first order difference is the quotient used to approximate the first derivate in the scent method iteration 3.6 MODELLING METHOD Step 1: obtain the langrange interpolation using

∏

Do this up to 5th order, the idea behind it, is to found a polynomial which agrees with special data. This polynomial can then be used to generate the number value at other pounds. We seek the polynomial of minimum degree which satisfies the interpolation condition P(

Step 2: use the polynomial to fixed the number values of the nodes Step 3: compare the found curves on the same graph and calculate the error bound used. Step 4: use finite difference, if it agrees Step 5: Transfer the graph to curve matching and see whether it is different/reasonable to use for interpretation.

3.7 FINITE DIFFERENCE TO CHECK THE DATA CONSISTENCY Forward difference: Shift operator: E =E

-

=

……………………….…………………..3.22

-

………………………………………….……………3.23

=

……………………………………………………………………3.24

Backward different: ∇ Central difference

= =

-

⁄

………………………….………..3.25

= -

⁄

=

⁄

…………………………….3.26

3.8 HIGHER ORDER DIFFERENCE We shall use the 1st order difference: equation 3.22 and 3.24 =

Generally;

-

=

=

…...3.27

-2

-

=∑

( )

…………………..….3.28

Generally: For a set of nodes i=0,1,…………………….N

CHAPTER FOUR 4.0 INTRODUCTION The study utilizes the efficiency of lagrangian interpolation and polyfit from MATLAB software. The curves obtained from the field and the generated one from Lagrange were compared to obtained good fitted curves for interpolation. The level of accuracy was determined using the numerical difference and central difference tables. The estimates were compared. 4.1 INTERPRETATION FROM THE LAGRANGE GENERATED CURVES The Lagrange in equation 3.1 was developed from 5th degree polynomial consisting of two variables in the form

( ) = ( ). This was used to estimate other points outside 1m,6m apart.

The curve generated when ran through a polyfit section of MATlabs programme,and also plotted by Techcal programme the curves showed an appreciable HK type from VES 38, in both cases when compared, the two Software showed an HK trend of two layer model. The error estimated is approximately at 3.3% and 2.5% respectively for MATlab and Techcal respectively.

Fig 4.0

Fig 4.1

VES 38 conventional curve for Electrical resistivity data.

VES 38 model curve for Electrical resistivity data generated by Matlab programme

Fig 4.2

VES 38 model curve for Electrical resistivity data generated by TechCal programme

The polynomial obtained within a, b  range is: -0.5083

+6.75

, which represents the trend of the whole data when compared with the whole values outside the range  a, b , used for the evaluation, the efficacy of the lagragian is that it reduces the task of plotting the whole data, the lagragian mimic the trend of the whole data obtained from the field, when this polynomial was interpreted on the master curves it gave about the same error values or even better in some of the VES points, which means that the lagragian has been able to filter out some of the noise signals obtained from the field due to certain factors. When the interpretation was carried out using the conventional methods in the three VES points, we recorded the error level of 2.5%, which is the same as one of the error obtained when lagragian is used.

Another curve was generated through section of MATlab, these points were plotted by the software and curve showed an appreciate KHK type for VES17 and also when plotted using Techcal programme. It also showed a HK three model curve. The error estimated at 6.9%.

Fig 4.3


Fig 4.4


Fig 4.5


The polyfit obtained is of the form: -4.38

+59.5

The above curves of fig4.12 and fig4.13 were transferred to the standard curves for segment by segment curve matching and resistivity interpreted by a log-linear graph. It gives the same number of layers with error level stated at 6.9% and 6.6% respectively.The error obtained from the conventional method on this VES point was 6.4%

In VES19 result was generated by repeating the process of 4.1 and 4.12 was not different as the error estimated is at 5.6% for Techcal and 12.9% from the MATlab respectively.

Fig 4.6


Fig 4.7


Fig 4.8


The polynomial obtained is -0.45

+6.5

-42.25

+156

-20.8 with HK type curves. When the

the conventional method was used, the error obtained was 5.8%, which is close to the polynomial obtained from the Techcal software. These two curves were transferred to a log-linear graph for interpretation it gives the same number of layers at two model layer with errors level stated at error 5.6% and respectively. 4.2 INTERPRETATION FROM FIELD DATA

The field data of the VES 38, when curve matched using win-Resist showed that the curve type is Hk type and has three layers with a bedrock extending to a value of about 245.5

With

maximum error of about 2.5% which is not far from what was obtained from the Lagrange equation 3.1, when the conventional method was used, but when the lagragian was used, it resulted in the same type of curve HK, with three layers an extending bedrock of 207.3 Also in VES 17 when the conventional method was used we have 4 layers, with a weathered basement of 12.6

When compared with the lagragian, using MATlab software we have 3

layers, with fresh basement of 98182.2

The result obtained when the Techcal was used gave

a 3 layer with weathered basement 18.5 In the VES point 19, using conventional method, three layers was obtained with a fresh basement of 1117.1

The result obtained from Techcal in very close to that of the conventional method

in the error values of 5.8% for conventional method and 5.6% for Techcal. The coventionalhas three layers, with high fresh bedrock value, but the Techcal has four layers with a wheared basement of 8.8

this last was not shown in the conventional method, however the MATlab

result showed an appreciable high level of error of about 12.9% with two layers an weathered bedrock of 82.9%.

4.3 DEDUCTIONS FROM THE LAGRANGE INTERPOLATION The lagrange showed that for a series of data obtain from any location, the lagrangian of a minimum order can be used to generate the curve which is a good representation of the whole data when compared with the field data and that p( ) =

= f( );

Such that )= obey

where the parameter =0i

,

= 1, i=j

If the data occurred within [a,b].

CHAPTER FIVE 5.0 CONCLUSION This research has examined the efficacy of Lagragian in the interpretation of electrical resistivity data. The results obtained are appreciate and close to the conventional system of interpretation, the generated curves from the Lagrangian interpolation showed a good trend which are representative of the intended patterns, even when the polynomial was generated with range [a,b] of the entire data. It shows that the Lagrangian when used to generated field curve at minimum degree, are intelligent to guest the trend of the whole data when considered as a whole at

[a,b]

generated polynomial. It’s also reduce the cost of time saving and curve matching stresses.

5.1 RECOMMENDATION In view of the above analysis, it is recommended that the project should be examined under different types of configurations such as: wenner, pole-dipole and double dipole arrays. Also, it is also recommended that the programmable aspect, which include the embedded system that can process and interpret the resistivity data as stated in the chapter one of this research.

REFERENCES American society of civil engineers, (1961), Ground-water basin management: Manual Eng. Practice. Archie, G.E, (1942), the electrical resistivity log as an aid in determining some reservoir characteristics: Petroleum Transactions of AIME. Bing, Z. and Greenhalgh, S.A., (2001), finite element three-dimensional direct current resistivity modeling Accuracy and efficiency considerations. Geophys.int. Brewitt-Taylor, C.R. and Weavear, J.T. (1976). On the finite difference solution of two dimensional induction problems. Geophys. J.R. astr. Soc. Broyden, C.G. (1972). Quasi-newton methods. In W.Murray, editor, Numerical methods for unconstrained optimization. Academics Press. Ellis, R.G. and Oldenburg,D.W. (1994). The pole-pole 3-dc resistivity inverse problem. A conjugate gradient approach Geophys.J.Int. Farquharson, C.G. and Oldenburg, D.W.(1998). Non-linear inversion using general measures of data misfit and model structure. Geophys.J.Int. Li, Y. and Oldenburg, D.W. (1992). Approximate inverse mappings in dc resistivity problems. Geophys. J. Int. Spitzer, K. (1998). Thre three-dimensional dc sensitivity for surface and subsurface sources. Geophys. J.Int Spitzer, K.(1999). Development of three-dimensional finite difference resistivity, sensitivity, and IP forward modeling techniques and their application to surface and borehole surveys. Habilitation thesis, University of Leipzig.

Zhang, J., Mackie, R.L., and Madden, T.R. (1995). 3-dc resistivity forward modeling and inversion using conjugate gradient. Geophys. J.Int. Zhao, S. and Yedlin, M.J. (1996). Some refinements on the finite difference method for 3-dc resistivity modeling. Geophys. J. Int.

APPENDIX

ELECTRICAL RESISTIVITY FIELD RECORD SHEET Project Name: NUMERICAL MODELING FOR ELECTRICAL RESISTIVITY DATA USING LAGRANGIAN INTERPOLATION Location: ALAGBAKA GRA

Instrument: TEREMETER RSO

Electrode Configuration: Schlumberger

Date: 29 – 01-2015

VES NO: VES 38 Co-ordinates: N0 station 1 2 3

AB/2 1 2 3

VES Description: 1 01.

E00 MN/2 0.25 0.25 0.25

1 49.

EL 356

K 6.28 25.12 56.52

R(ohm) 12.3636 2.9650 1.1678

App.R 78 74 66

V 1768 424 167

I 143 143 143

4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

4 6 6 8 12 15 15 25 32 40 40 65 100 100 150 225 225 325 450

0.25 0.25 0.50 0.50 0.50 0.50 1.00 1.00 1.00 1.00 2.50 2.50 2.50 5.00 5.00 5.00 10.00 10.00 10.00

100.48 226.08 113.04 200.96 452.16 706.50 353.25 981.25 1607.68 2512.00 1004.80 2653.30 6280.00 3140.00 7065.00 15896.25 7948.13 16583.13 31792.50

0.6408 0.2535 0.6690 0.3472 0.1389 0.09412 0.1957 0.0956 0.0672 0.0556 0.1181 0.06812

64 57 76 70 63 67 69 94 108 140 119 181

91 36 95 50 20 13 27 13 9 8 17 9

142 142 142 144 144 138 138 136 134 144 144 132




Date: 29 – 01-2015

VES NO: VES 38 Co-ordinates: N0 station 1 2

AB/2 1 2

VES Description: 1 00.

E00 MN/2 0.25 0.25

1 17. K 6.28 25.12

EL 332 R(ohm) 15.6867 6.6304

App.R 99 167

V 2353 1915

I 150 138

3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

3 4 6 6 8 12 15 15 25 32 40 40 65 100 100 150 225 225 325 450

0.25 0.25 0.25 0.50 0.50 0.50 0.50 1.00 1.00 1.00 1.00 2.50 2.50 2.50 5.00 5.00 5.00 10.00 10.00 10.00

56.52 100.48 226.08 113.04 200.96 452.16 706.50 353.25 981.25 1607.68 2512.00 1004.80 2653.30 6280.00 3140.00 7065.00 15896.25 7948.13 16583.13 31792.50

3.6463 2.2662 0.9493 2.2464 1.2788 0.4058 0.22 0.4133 0.14 0.1095 0.07692 0.1329 0.1515 0.04918 0.1

206 228 215 254 257 183 155 146 137 176 193 134 402 309 314

536 315 131 310 113 56 33 62 21 15 11 19 20 6 12

147 139 138 138 104 138 150 150 150 137 143 143 132 122 120




Date: 29 – 01-2015

VES NO: VES 17

VES Description:

Co-ordinates: N0 station 1

AB/2 1

1 52.

E00 MN/2 0.25

1 34. K 6.28

EL 340.2 R(ohm) 20.7286

App.R 130

V 29902

I 140

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

2 3 4 6 6 8 12 15 15 25 32 40 40 65 100 100 150 225 225 325 450

0.25 0.25 0.25 0.25 0.50 0.50 0.50 0.50 1.00 1.00 1.00 1.00 2.50 2.50 2.50 5.00 5.00 5.00 10.00 10.00 10.00

25.12 56.52 100.48 226.08 113.04 200.96 452.16 706.50 353.25 981.25 1607.68 2512.00 1004.80 2653.30 6280.00 3140.00 7065.00 15896.25 7948.13 16583.13 31792.50

6.027 2.146 1.3301 0.6136 1.431 0.6148 0.2786 0.2049 0.4083 0.2258 0.1719 0.1071 0.3841 0.0753 0.0154 0.0231

151 121 134 139 162 124 126 145 144 222 276 269 386 200 97 72

675 279 173 184 189 83 39 25 49 28 22 15 53 11 2 3

MATlab

VES 38 h1=12.7 h2=13.5 h3=21.4 h4=54.3

VES 17

Error = 3.3% HK CURVES

112 130 130 132 132 135 140 122 120 124 128 140 138 146 130 130

h1=9.3 h2=11.9 h3=26.6

HKH CURVES

Error = 6.9%

K CURVES

Error = 12.9%

VES19 h1=0.4 h2=246.7

Techcal Software

VES 38 h1=0.4 h2=14.7 h3=27.2


VES 17 h1=0.6 h2=3.8 h3=5.3 h4=120.

KHK CURVES

Error = 5.6%

K CURVES

Error = 12.9%

VES19 h1=0.4 h2=246.7

Conventional

VES 38 h1=0.9 h2=3.8 h3=9.4


VES 17 h1=0.9 h2=1.6 h3=6.6

KHK CURVES

Error = 6.4%

KH CURVES

Error = 5.8%

VES19 h1=0.6 h2=2.

MATLAB SCRIPT %progam to simulate polynomial lagrangian by aminu I.oluwagbenga x=[ 1 2 3 4 6] y=[99 167 206 228 215] %solve lagrangian polynomial [P,R,S]=lagrangepoly(x,y) xx=0.5:0.01:6.5 %plot the polynomial fitted curve plot(xx,polyval(p,xx),x,y) grid axis([0.5 8.5 -5 5])

VES 19

VES 38

VES 17

NUMERICAL MODELING OF ELECTRICAL RESISTIVITY DATA USING LAGRANGIAN INTERPOLATION

Recommend Documents