How to generate a diffusing resistive structure to do forward modelling?

[jcmefra]

Problem description

I don't know if it's possible to create a mesh with a diffusing resistive structure, I want to model an aquifer with geothermal water infiltration so the mesh would be basically one layer (aquifer) + diffusing structure inside (geothermal water infiltration)

Your environment

---

Date: Fri Aug 25 23:54:59 2023 -05

            OS : Linux
        CPU(s) : 8
       Machine : x86_64
  Architecture : 64bit
           RAM : 15.6 GiB
   Environment : Jupyter
   File system : ext4

Python 3.10.12 | packaged by conda-forge | (main, Jun 23 2023, 22:40:32) [GCC 12.3.0]

       pygimli : 1.4.3
        pgcore : 1.4.0
         numpy : 1.25.0
    matplotlib : 3.7.2
         scipy : 1.11.2
       IPython : 8.14.0
       pyvista : 0.41.1

---

Steps to reproduce

It's more like a question, not an error.

Expected behavior

To do a mesh with a diffusing structure

Actual behavior

I don't know how to do it

[carsten-forty2]

I don't get it completely. Should be the mesh structure (geometry) itself time depending? Or do you want field values for the mesh time depending? Do you have information about the infiltration, e.g. like from a fluid flow simulator?

[jcmefra]

I'm not that good at PyGimli since I'm just learning how to use it, I was checking the documentation and it shows how to do forward modelling for 2D ERT, however it only shows how to do a basic modelling (shape, markers, rho value for each marker) and also how to add topography data. I want to know if it's possible to use pygimli to create a diffusing shape and assign it a marker too. Thanks!

Edit: I don't want the structure time depending, just to have diffusive borders as it's water infiltration so it's not a defined body.

[halbmy]

I guess we still don't understand what you mean with diffusive structure. Do you mean a statistical behaviour like here? https://www.pygimli.org/_tutorials_auto/3_inversion/plot_6-geostatConstraints.html#generating-geostatistical-media

[jcmefra]

I guess we still don't understand what you mean with diffusive structure. Do you mean a statistical behaviour like here? https://www.pygimli.org/_tutorials_auto/3_inversion/plot_6-geostatConstraints.html#generating-geostatistical-media

Hello, thanks for the information. Yes it can be something like that, or maybe something like hydrogeophysical modelling but it's more complex and I don't want it to be time depending, just two dimensional (x, depth).

Here's my plot:

And here are the resistivities I defined:

It can work(ish) but I want something more precise por that shape defined by 4, 5 and 6 markers as it's supposed to be a fluid infiltrating into the first layer.

[carsten-forty2]

For pure diffusive behaviour you can try the script below. It calculates the diffusion of a normalized concentration from an injection point. There are the stationary concentration at infty times (see screenshot) and also some evolution over time so you can choose a time step that fits your need (animation at the end). The boundary conditions are kind of arbitrary, just to show the way. In this case, there is outflow of concentration in the subsurface and fixed concentration of zero on the surface. (Edited to more resonable results) You can switch it below but the result may be a bit more boring, However, its fun to play arround a little.

The concentration can be translated into a resisitvity (or better a conductivity) offset for ERT simulation as ever you like.

If you want the mesh structure, i.e. the cell edges, follow the contour lines would be complicated but not impossible. But I doubt its needed for an ERT simulation.

If you want a more plume like structure, then you need an additional advective term which is currently not yet possible with the current gimli state but work in progress. However, you can try to fake it a little more with some more boundary condition tweaking.

import numpy as np
import pygimli as pg

import pygimli.meshtools as mt

world = mt.createWorld(start=(-350, -250), end=(350, 0))

iPos = (-120, -200)
layer = mt.createPolygon([(-350, -150), iPos, (350, -150)], 
                        addNodes=100, interpolate='spline', boundaryMarker=5)

cut = mt.createPolygon([iPos, (-180, -250)])

world += layer
world += cut
world.addRegionMarker((300, -10), marker=2, area=100)

pg.show(world, markers=True, showNodes=True)

mesh = mt.createMesh(world)
#pg.show(mesh, markers=True, showMesh=True)

diff = pg.solver.cellValues(mesh, {0: 1e-2, 1: 1e-2, 2: 100})
pg.show(mesh, diff, label='diffusity')

iPosID = mesh.findNearestNode(iPos)

# stationary solution 
C = pg.solver.solve(mesh, a=diff,
                    bc={'Dirichlet': {'-2': 0.0}, 'Neumann': {'-1': -1},
                        'Node':[iPosID, 100]}, verbose=True)

pg.show(mesh, C, label='concentration', showMesh=True, 
        cMin=0, cMax=100,  nCols=10, nLevs=11, linewidths=0.5)

# evolution over time
C = pg.solver.solve(mesh, a=diff,
                    bc={'Dirichlet': {'-1': 0.0}, 'Neumann': {'-2': -1},
                        'Node':[iPosID, 100]}, 
                    times=np.linspace(0,100,100), verbose=True)

# show animation of every 2nd timestep
pg.show(mesh, C[::2], label='concentration', showMesh=True, 
        cMin=1, cMax=100,  nCols=10, nLevs=11, linewidths=0.5)

Edit 2.: Little bit more plume like with anisotrop diffusity:

diff = pg.solver.cellValues(mesh, {0: pg.solver.createAnisotropyMatrix(0.1, 0.1, 0.0),
                                   1: pg.solver.createAnisotropyMatrix(0.1, 0.1, 0.0),
                                   2: pg.solver.createAnisotropyMatrix(1, 100, 20.0*np.pi/180)}
                                   )

Time depending simulation seems broken with anisotrope material but I guess steady simulation should be sufficient here.

[jcmefra]

For pure diffusive behaviour you can try the script below. It calculates the diffusion of a normalized concentration from an injection point. There are the stationary concentration at infty times (see screenshot) and also some evolution over time so you can choose a time step that fits your need (animation at the end). The boundary conditions are kind of arbitrary, just to show the way. In this case, there is outflow of concentration in the subsurface and fixed concentration of zero on the surface. (Edited to more resonable results) You can switch it below but the result may be a bit more boring, However, its fun to play arround a little.

The concentration can be translated into a resisitvity (or better a conductivity) offset for ERT simulation as ever you like.

If you want the mesh structure, i.e. the cell edges, follow the contour lines would be complicated but not impossible. But I doubt its needed for an ERT simulation.

If you want a more plume like structure, then you need an additional advective term which is currently not yet possible with the current gimli state but work in progress. However, you can try to fake it a little more with some more boundary condition tweaking.

import numpy as np
import pygimli as pg

import pygimli.meshtools as mt

world = mt.createWorld(start=(-350, -250), end=(350, 0))

iPos = (-120, -200)
layer = mt.createPolygon([(-350, -150), iPos, (350, -150)], 
                        addNodes=100, interpolate='spline', boundaryMarker=5)

cut = mt.createPolygon([iPos, (-180, -250)])

world += layer
world += cut
world.addRegionMarker((300, -10), marker=2, area=100)

pg.show(world, markers=True, showNodes=True)

mesh = mt.createMesh(world)
#pg.show(mesh, markers=True, showMesh=True)

diff = pg.solver.cellValues(mesh, {0: 1e-2, 1: 1e-2, 2: 100})
pg.show(mesh, diff, label='diffusity')

iPosID = mesh.findNearestNode(iPos)

# stationary solution 
C = pg.solver.solve(mesh, a=diff,
                    bc={'Dirichlet': {'-2': 0.0}, 'Neumann': {'-1': -1},
                        'Node':[iPosID, 100]}, verbose=True)

pg.show(mesh, C, label='concentration', showMesh=True, 
        cMin=0, cMax=100,  nCols=10, nLevs=11, linewidths=0.5)

# evolution over time
C = pg.solver.solve(mesh, a=diff,
                    bc={'Dirichlet': {'-1': 0.0}, 'Neumann': {'-2': -1},
                        'Node':[iPosID, 100]}, 
                    times=np.linspace(0,100,100), verbose=True)

# show animation of every 2nd timestep
pg.show(mesh, C[::2], label='concentration', showMesh=True, 
        cMin=1, cMax=100,  nCols=10, nLevs=11, linewidths=0.5)

Edit 2.: Little bit more plume like with anisotrop diffusity:

diff = pg.solver.cellValues(mesh, {0: pg.solver.createAnisotropyMatrix(0.1, 0.1, 0.0),
                                   1: pg.solver.createAnisotropyMatrix(0.1, 0.1, 0.0),
                                   2: pg.solver.createAnisotropyMatrix(1, 100, 20.0*np.pi/180)}
                                   )

Time depending simulation seems broken with anisotrope material but I guess steady simulation should be sufficient here.

May I know how to do the inversion? I replaced Conc for Rho and did some workaround to have iso-resistivity curves:

rho = pg.solver.solveFiniteElements(mesh, a=diff,
                    bc={'Dirichlet': {'-2': 100}, 'Neumann': {'-1': 1},
                        'Node':[iPosID, 1]}, verbose=True)

pg.show(mesh, rho, label='Resistivity', showMesh=True, 
        cMin=0, cMax=100,  nCols=10, nLevs=11, linewidths=0.5)

I did this:

scheme = ert.createData(elecs=np.linspace(start=-350, stop=350, num=96),
                           schemeName='slm')

data = ert.simulate(mesh, scheme=scheme, res=rho, noiseLevel=3,
                    noiseAbs=1e-6, seed=1337)

pg.info(np.linalg.norm(data['err']), np.linalg.norm(data['rhoa']))
pg.info('Simulated data', data)
pg.info('The data contains:', data.dataMap().keys())

pg.info('Simulated rhoa (min/max)', min(data['rhoa']), max(data['rhoa']))
pg.info('Selected data noise %(min/max)', min(data['err'])*100, max(data['err'])*100)

But it says "Simulate called with wrong resistivity array.", because res needs to be in format [[marker, rho], [marker2, rho2],...]. Is there a way to convert that into polygons with markers, maybe?

Thanks in advance for your help!

[carsten-forty2]

The unit of your rho is no resistivity but a normalized concentration .. you need to translate it somehow. See example below.

For ERT simulation you need a cell based value array so you have to interpolate the node based values into cell based values. It should be possible to pass this array into ert.simulate directly. (In some further versions of pg such simulations can be done with nodes based values directly.)

Note, the linear relation from any concentration to resistivity given below is not related to any physical behaviour. You should do some paper search to find a better function.

For meaningfull accuracy of your ERT simulation you also need to add some larger boundary as well as some local refinement at the electrode positions. Please consult the appropriate ert examples and manuals.

# Background rho0
rhoMin = 2
rhoMax = 200
rho0 = pg.solver.cellValues(mesh, {0: 500,
                                   1: rhoMin,
                                   2: rhoMax}
                                   )

pg.show(mesh, rho0, label='rho background')
C = pg.interpolate(mesh, C, mesh.cellCenters())
C[C < 0] = 0

# add anomal resistivity as linear function from rhoMin to rhoMax depending on concentration
rho = 1/(1/np.array(rho0) + 1/(rhoMin)*(C/100))

pg.show(mesh, rho, label='rho', logScale=True)

[jcmefra]

The unit of your rho is no resistivity but a normalized concentration .. you need to translate it somehow. See example below.

For ERT simulation you need a cell based value array so you have to interpolate the node based values into cell based values. It should be possible to pass this array into ert.simulate directly. (In some further versions of pg such simulations can be done with nodes based values directly.)

Note, the linear relation from any concentration to resistivity given below is not related to any physical behaviour. You should do some paper search to find a better function.

For meaningfull accuracy of your ERT simulation you also need to add some larger boundary as well as some local refinement at the electrode positions. Please consult the appropriate ert examples and manuals.

# Background rho0
rhoMin = 2
rhoMax = 200
rho0 = pg.solver.cellValues(mesh, {0: 500,
                                   1: rhoMin,
                                   2: rhoMax}
                                   )

pg.show(mesh, rho0, label='rho background')
C = pg.interpolate(mesh, C, mesh.cellCenters())
C[C < 0] = 0

# add anomal resistivity as linear function from rhoMin to rhoMax depending on concentration
rho = 1/(1/np.array(rho0) + 1/(rhoMin)*(C/100))

pg.show(mesh, rho, label='rho', logScale=True)

Got it, thanks! One las thing, rhomap for ERT simulation must have pairs of [marker, rhovalue], but I only have 3 region markers and the decreasing rho values should be into the first layer/marker. What should I do?

Another little question: I'm trying to do the ERT inversion with the world geometry and the mesh you provided above (ignoring the diffusive fluid as I don't know how to add it yet, just 3 shapes with 3 markers and 3 rho values), I'm getting a very poor pseudosection, so maybe I'm not going to be able to run the inversion properly. How to address that? I defined my markers to 1, 2, 3:

world = mt.createWorld(start=(-350, -250), end=(350, 0))
iPos = (-125, -200) #Inyection Point

fault = mt.createPolygon([iPos, (-150, -250)], isClosed= True)

layer = mt.createPolygon([(-350, -150), iPos, (350, -150)], 
                         addNodes=100, interpolate='spline', isClosed=False)

world += layer
world += fault

world.addRegionMarker((300, -10), marker=1)
world.addRegionMarker((-300, -200), marker=2)
world.addRegionMarker((300, -200), marker=3)

pg.show(world)

scheme = ert.createData(elecs=np.linspace(start=-350, stop=350, num=96),
                           schemeName='slm')

mesh = mt.createMesh(world, quality = 34)

for p in scheme.sensors():
    world.createNode(p)
    world.createNode(p - [0, 0.1])

rhomap = [[1, 80.],
          [2, 100.],
          [3, 500.]]

pg.show(mesh, data=rhomap, label=pg.unit('res'), showMesh=True)

data = ert.simulate(mesh, scheme=scheme, res=rhomap, noiseLevel=3,
                    noiseAbs=1e-6, seed=1337)

pg.info(np.linalg.norm(data['err']), np.linalg.norm(data['rhoa']))
pg.info('Simulated data', data)
pg.info('The data contains:', data.dataMap().keys())

pg.info('Simulated rhoa (min/max)', min(data['rhoa']), max(data['rhoa']))
pg.info('Selected data noise %(min/max)', min(data['err'])*100, max(data['err'])*100)

data.remove(data['rhoa'] < 0)
pg.info('Filtered rhoa (min/max)', min(data['rhoa']), max(data['rhoa']))

data.save('Forward_mod.dat')

ert.show(data)

Thanks!

Edit: when I define the geometry as world + polygons using poly = mt.createPolygon(..., marker = n, markerPosition = coords) and then combining them by doing geom = world + poly1 + poly2 + poly3 + ... it works, and I'm able to run the inversion. I defined world and 5 different polygons that way and combined them to do this:

And my inversion was okay:

However, that way the boundary markers are pretty broken so I can't simulate the diffusive fluid as you did:

Note: The polygons 4, 5 and 6 (the "plume") are going to be replaced by the diffusing body.

Thanks in advance.

[carsten-forty2]

Are there any questions remain? Can't see them in this multi edited thread ;)

Just one additional remark, you can pass the rho-per-cell-data array directly into the ert simulation. A region to value map is only the default if you don't have the complete vector.

[carsten-forty2]

ah .. regarding the last image .. the outside boundary marker are important for ERT simulation. -1 is allways used as surface and -2 is allways used as subsurface. Values > 0 are ignored. In your case with the last image .. the lower boundaries with marker == 1 are wrong for ERT simulation.

[jcmefra]

Are there any questions remain? Can't see them in this multi edited thread ;)

Just one additional remark, you can pass the rho-per-cell-data array directly into the ert simulation. A region to value map is only the default if you don't have the complete vector.

Just a couple:

1. why my pseudo section (ignoring the plume, just simulating the ERT data based on the layers) is too poor when I define the geometry like you did? I mean world + layer + cut + add region markers. If I define my geometry like that the boundaries are as expected, but somehow the ERT simulated data is very poor. Am I missing something?

2. I try to add "rho" vector (calculated based on C) to my ert simulate but it says the res data is wrong. I like how my plume looks, now I need to fix the ERT data simulation issues. Thanks!

[carsten-forty2]

Okay .. you issueing two problems.

1. missing appropriate local refinement at the electrode positions. The mesh need to be created after you add them to the world geometry. Its good to have at least on cell for each sensor.
2. missing appropriate background region for forward modelling

The mesh we created in the first should then be used for the ert simulation. For the inversion you need then the other mesh to avoid an inverse crime.

See the attached script addressing these both issues without further optimization. You should play with the parameters to become familiar with the effects.

import numpy as np
import pygimli as pg

import pygimli.meshtools as mt
from pygimli.physics import ert

world = mt.createWorld(start=(-350, -250), end=(350, 0))

iPos = (-120, -200)
layer = mt.createPolygon([(-350, -150), iPos, (350, -150)], 
                        addNodes=100, interpolate='spline', boundaryMarker=5)

cut = mt.createPolygon([iPos, (-180, -250)])

world += layer
world += cut

world.addRegionMarker((300, -10), marker=2, area=500)

# Sensors should no be on the corner so we put them a little inside
scheme = ert.createData(elecs=np.linspace(start=-330, stop=330, num=96),
                           schemeName='slm')

# we need local refinement at the electrodes to achieve sufficient accuracy
for p in scheme.sensors():
    world.createNode(p)
    world.createNode(p - [0, 0.1])

pg.show(world, markers=True, showNodes=True)

mesh = mt.createMesh(world, quality=33, area=100)

# we need larger boundary to avoid numerical issues
mesh = mt.appendTriangleBoundary(mesh, xbound=1000, ybound=500, marker=3)
# Note. appendTriangleBoundary is just an easy fallback. You will get better results if you create the boundary manually on geometry level.

pg.show(mesh, markers=True, showMesh=True)

rhoMin = 2
rhoMax = 200
rho0 = pg.solver.cellValues(mesh, {0: 500,
                                   1: rhoMin,
                                   2: rhoMax,
                                   3: 100,}
                                   )

pg.show(mesh, rho0, label='rho background', showMesh=True)

# we simulate the ert with our forward mesh
data = ert.simulate(mesh, scheme=scheme, res=rho0, 
                    noiseLevel=3, noiseAbs=1e-6, seed=1337)

pg.show(data)

diff = pg.solver.cellValues(mesh, {'0,1,3': pg.solver.createAnisotropyMatrix(0.1, 0.1, 0.0),
                                   2: pg.solver.createAnisotropyMatrix(1, 100, 20.0*np.pi/180)}
                                   )

iPosID = mesh.findNearestNode(iPos)

# stationary solution 
C = pg.solver.solve(mesh, a=diff,
                    bc={'Dirichlet': {'-1': 0.0}, 'Neumann': {'-2': -1},
                        'Node':[iPosID, 100]}, verbose=True)

pg.show(mesh, C, label='concentration', showMesh=True, 
        cMin=0, cMax=100, nCols=10, nLevs=11, linewidths=0.5)

C = pg.interpolate(mesh, C, mesh.cellCenters())
C[C < 0] = 0

# add anomal resistivity as linear function from rhoMin to rhoMax depending on concentration
rho = 1/(1/np.array(rho0) + 1/(rhoMin)*(C/100))

data = ert.simulate(mesh, scheme=scheme, res=rho, 
                    noiseLevel=3, noiseAbs=1e-6, seed=1337)

pg.show(data)

[jcmefra]

diff = pg.solver.cellValues(mesh, {0: 1e-2, 1: 1e-2, 2: 100})

So I was having issues with the mesh. Thank you so much for your help! I really appreciate it. Have a great day.