Question about petrophysical individually inversion.

[shdlovezxm]

Problem description

hello, "When I assume ertTrans=ArchieTrans (rFluid=20, phi=0.3), can I directly obtain the Saturation model from ERT data? The following is my code, and I did not get the desired result. "

Please help me, thank you very much.

Steps to reproduce

Tell us how to reproduce this issue. Ideally, you could paste the code that produces the error:

import matplotlib.pyplot as plt

import pygimli as pg
from pygimli.physics import ert
import pygimli.meshtools as mt
from pygimli.physics.petro import transFwdArchieS as ArchieTrans
from pygimli.frameworks import PetroInversionManager

data = ert.load("./lake.ohm")
print(data)

data['k'] = ert.createGeometricFactors(data, numerical=True)
print(data)

k0 = ert.createGeometricFactors(data)

ert.showData(data, vals=k0/data['k'], label='Topography effect')

mgr = ert.ERTManager(data)

mgr.checkData(data)
pg.show(data)

print(data)

data['err'] = ert.estimateError(data, relativeError=0.02, absoluteUError=100e-6)
print(data)
# or manually:
# data['err'] = data_errors  # somehow
ert.show(data, data['err']*100, label="error (%)")

inv = mgr.invert(data, lam=20, verbose=True, quality=33.6)
mgr.showResult()

pg.info("ERT Petrogeophysical Inversion")
ertTrans = ArchieTrans(rFluid=20, phi=0.3)
ERTPetro = PetroInversionManager(petro=ertTrans, mgr=mgr)
satERT = ERTPetro.invert(data, limits=[0., 1.], lam=20,
                         verbose=True)
ERTPetro.inv.echoStatus()
satKW = dict(logScale=False, cMap="plasma_r")
pg.show(ERTPetro.paraDomain, ERTPetro.paraModel(satERT), **satKW, 
             showMesh=True, label=r'Saturation (${\tt petro}$)')

...

Actual behavior

Tell us what happens instead and/or provide the output of your script.

23/10/23 - 15:34:31 - pyGIMLi - INFO - ERT Petrogeophysical Inversion
23/10/23 - 15:34:31 - pyGIMLi - INFO - Found 1 regions.
23/10/23 - 15:34:31 - pyGIMLi - INFO - Found 1 regions.
23/10/23 - 15:34:31 - pyGIMLi - INFO - Creating forward mesh from region infos.
23/10/23 - 15:34:31 - pyGIMLi - INFO - Creating refined mesh (H2) to solve forward task.
23/10/23 - 15:34:31 - pyGIMLi - INFO - Mesh for forward task: Mesh: Nodes: 1505 Cells: 2768 Boundaries: 2196
23/10/23 - 15:34:31 - pyGIMLi - INFO - Use median(data values)=42.69938678956654
23/10/23 - 15:34:31 - pyGIMLi - INFO - Created startmodel from forward operator: 692 [2.2813027502697536,...,2.2813027502697536]
23/10/23 - 15:34:31 - pyGIMLi - INFO - Starting inversion.
fop: <pygimli.frameworks.modelling.PetroModelling object at 0x000002DC0A0C8E00>
Data transformation: <pgcore._pygimli_.RTrans object at 0x000002DC0A1893A0>
Model transformation (cumulative):
     0 <pgcore._pygimli_.RTransLogLU object at 0x000002DC0A0A3A00>
min/max (data): 11.44/87.51
min/max (error): 2%/2.58%
min/max (start model): 2.28/2.28
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
inv.iter 1 ... chi² = 99372.5 (dPhi = -7909.13%) lam: 20
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 99372.5 (dPhi = 0.0%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 99372.5 (dPhi = 0.0%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 99372.5 (dPhi = 0.0%) lam: 20.0
################################################################################
#                 Abort criteria reached: dPhi = 0.0 (< 2.0%)                  #
################################################################################

[halbmy]

Your startmodel is somehow automatically determined with 2.28, which is outside of your given saturation limit [0, 1]. Can you try with a valid starting model invert(limits=[0, 1], startModel=0.5, ...)?

[shdlovezxm]

lake.txt

The range of the Saturation model is [0,1]. No matter how I change the startModel, I cannot obtain the desired result. If I remove the model range, I can obtain the result. However, the range of the Saturation model is [0.81, 5.35], which is the wrong result?

Actual behavior

satERT = ERTPetro.invert(data, lam=20, limits=[0, 1], startModel=0.5, verbose=True)

fop: <pygimli.frameworks.modelling.PetroModelling object at 0x000002DC84675EA0>
Data transformation: <pgcore._pygimli_.RTrans object at 0x000002DC8465FD60>
Model transformation (cumulative):
     0 <pgcore._pygimli_.RTransLogLU object at 0x000002DC8464D340>
min/max (data): 11.44/87.51
min/max (error): 2%/2.58%
min/max (start model): 0.5/0.5
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
inv.iter 1 ... chi² = 402481.45 (dPhi = 78.89%) lam: 20
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 132241.06 (dPhi = 67.14%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 99527.0 (dPhi = 24.7%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 99372.5 (dPhi = 0.22%) lam: 20.0
################################################################################
#                 Abort criteria reached: dPhi = 0.22 (< 2.0%)                 #
################################################################################

satERT = ERTPetro.invert(data, lam=20, verbose=True)

fop: <pygimli.frameworks.modelling.PetroModelling object at 0x000002DC9A609900>
Data transformation: <pgcore._pygimli_.RTrans object at 0x000002DC9A5D6D00>
Model transformation: <pgcore._pygimli_.RTransLog object at 0x000002DC9A609A90>
min/max (data): 11.44/87.51
min/max (error): 2%/2.58%
min/max (start model): 2.28/2.28
--------------------------------------------------------------------------------
--------------------------------------------------------------------------------
inv.iter 1 ... chi² = 94.54 (dPhi = 92.29%) lam: 20
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 11.3 (dPhi = 87.39%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 0.6 (dPhi = 87.54%) lam: 20.0

################################################################################
#                  Abort criterion reached: chi² <= 1 (0.60)                   #
################################################################################

[halbmy]

You took the lake case (https://www.pygimli.org/_tutorials_auto/3_inversion/plot_8-regionWise.html) which is an underwater case but you don't treat it like that. In the pure resistivity inversion, you obtain subsurface resistivities of about 8 Ohmm, which is 2.5 times below the specified fluid resistivity of 20 Ohmm. If you formally convert 8 Ohmm into saturations you end up with values above 1. So technically all is working well but the data set does not make sense in terms of Archie unless you include the water body. (In this case there is a lot of clay and Archie won't work anyway)

[shdlovezxm]

I understand. Thank you for your patient explanation.

[halbmy]

I guess it can be closed. If there is still a problem, just reopen.