Time-lapse inversion

In this tutorial, we will see how to use ResIPy API to time-lapse inversion of the problem in Figure 5.21 Note that more datasets are included here (as per Musgrave and Binley, 2011, doi:10.1016/j.jhydrol.2010.11.008

In [1]:
from resipy import R2
k = R2(typ='R2')
import matplotlib.pyplot as plt  # used for more control of plots
C:\Program Files (x86)\WPy64-3770\python-3.7.7.amd64\lib\site-packages\resipy\meshTools.py:50: UserWarning: pyvista not installed, 3D meshing viewing options will be limited
  warnings.warn('pyvista not installed, 3D meshing viewing options will be limited')
API path =  C:\Program Files (x86)\WPy64-3770\python-3.7.7.amd64\lib\site-packages\resipy
ResIPy version =  2.2.2
cR2.exe found and up to date.
R3t.exe found and up to date.
cR3t.exe found and up to date.
Working directory is: C:\Program Files (x86)\WPy64-3770\python-3.7.7.amd64\lib\site-packages\resipy
clearing dirname

Time-lapse datasets are stored in a sub-folder

In [2]:
k.createTimeLapseSurvey(r'C:\ResIPy\Notebook examples\Time-lapse inversion\data', ftype='ProtocolDC')
14/14 imported
In [3]:
k.importElec(r'C:\ResIPy\Notebook examples\Time-lapse inversion\electrodes.csv') # read electrode geometry
In [4]:
k.createMesh(typ='trian', cl=0.25, fmd=5) # triangular mesh with element size near electrodes equal to 0.25m. fmd sets depth for fine mesh
k.zlim = [-5, 0]
k.showMesh(xlim=[0,63],zlim=[-5,0])
Creating triangular mesh...Reading mesh.msh
Gmsh version == 3.x
reading node coordinates...
Determining element type...Triangle
Reading connection matrix...
ignoring 0 elements in the mesh file, as they are not required for R2/R3t
Finished reading .msh file
ResIPy Estimated RAM usage = 0.067461 Gb
done
In [5]:
k.param['b_wgt'] = 0.02 # we are now informing R2 that the relative error is 0.02 for inversion
k.param['a_wgt'] = 0.001 # the offset error is 0.001 ohms
In [6]:
k.invert(parallel=True)  # run time-lapse difference inversion in parallel once reference dataset is inverted
Writing .in file and protocol.dat... Matching quadrupoles between surveys for difference inversion...397 in common...done in 0.045907s
done!
------------ INVERTING REFERENCE SURVEY ---------------


 >> R  2    R e s i s t i v i t y   I n v e r s i o n   v4.0 <<

 >> D a t e : 11 - 08 - 2020
 >> My beautiful survey                                                             
 >> I n v e r s e   S o l u t i o n   S e l e c t e d <<
 >> Determining storage needed for finite element conductance matrix
 >> Generating index array for finite element conductance matrix
 >> Reading start resistivity from res0.dat                      
 >> R e g u l a r i s e d   T y p e <<
 >>   L i n e a r    F i l t e r    <<
 >> L o g - D a t a   I n v e r s i o n <<
 >> N o r m a l   R e g u l a r i s a t i o n <<
 >> D a t a   w e i g h t s   w i l l   b e  m o d i f i e d <<


 Processing dataset   1


 Measurements read:   397     Measurements rejected:     0
   Geometric mean of apparent resistivities:  0.78168E+02

 >> Total Memory required is:          1.958 Gb

   Iteration   1
     Initial RMS Misfit:        17.21       Number of data ignored:     0
     Alpha:          86.273   RMS Misfit:        3.12  Roughness:       28.655
     Alpha:          40.044   RMS Misfit:        3.10  Roughness:       40.278
     Alpha:          18.587   RMS Misfit:        3.21  Roughness:       55.497
     Step length set to   1.00000
     Final RMS Misfit:        3.10
     Updated data weights

   Iteration   2
     Initial RMS Misfit:         2.60       Number of data ignored:     0
     Alpha:          23.933   RMS Misfit:        1.30  Roughness:       47.766
     Alpha:          11.109   RMS Misfit:        0.94  Roughness:       68.088
     Step length set to   1.00000
     Final RMS Misfit:        0.94
     Final RMS Misfit:        1.00

 Solution converged - Outputing results to file

 Calculating sensitivity map


 Processing dataset   2


 End of data:  Terminating
--------------------- MAIN INVERSION ------------------
________________System-Check__________________
Kernel type: Windows
Processor info: Intel64 Family 6 Model 142 Stepping 9, GenuineIntel
4 Threads at <= 2401.0 Mhz
Total memory = 15.7 Gb (usage = 47.0)
Wine Version = Native Windows (N/A)
13/13 inversions completed
----------- END OF INVERSION IN // ----------
14/14 results parsed (14 ok; 0 failed)
In [7]:
k.meshResults[0].df.columns # this will display the attributes that can be plotted
Out[7]:
Index(['param', 'elm_id', 'region', 'cellType', 'X', 'Y', 'Z',
       'Resistivity(ohm.m)', 'Resistivity(log10)', 'Sensitivity(log10)',
       'Conductivity(mS/m)'],
      dtype='object')

Plot log resistivity image for first (reference) dataset and show sensitivity pattern (see Figure 5.23b)

In [8]:
fig, axs = plt.subplots(2, 1, figsize=(18,6))
ax=axs[0]
k.showResults(ax=ax, index=0, attr='Resistivity(log10)', contour=False, vmin=1.5, vmax=2.4, color_map='jet', cropMaxDepth=True)
ax.set_title('09 Feb 2005  log10 resistivity')
# We can now plot the sensitivity map (as shown in Figure 5.23b)
ax=axs[1]
k.showResults(ax=ax, index=0, attr='Sensitivity(log10)', contour=True, vmin=-3, vmax=0, color_map='jet', cropMaxDepth=True)
ax.set_title('Sensitivity plot')
Out[8]:
Text(0.5, 1.0, 'Sensitivity plot')

Now plot resistivity difference for the subsequent datasets relative to first dataset

In [10]:
headings = ['10 Mar 2005','07 Apr 2005','21 Apr 2005','03 May 2005','17 May 2005', '07 Jun 2005', '28 Jun 2005', '12 Jul 2005', '30 Aug 2005', '20 Sep 2005','30 Nov 2005', '18 Jan 2006', '09 Feb 2006']
fig, axs = plt.subplots(13, 1, figsize=(18,30))
for i in range(len(headings)):
    ax=axs[i]
    k.showResults(ax=axs[i], index=i+1, attr='difference(percent)', contour=False, vmin=-30, vmax=30, color_map='seismic', cropMaxDepth=True)
    ax.set_title(headings[i])
    if i < 12:
        ax.set_xlabel('')
In [ ]: