Extraction well simulation

[mariosgeo]

Problem description

This is not a code issue, rather seeking advice on the modeling. I am trying to simulate a simple extraction well, for various scenarios (no ground water flow, some ground water flow etc). The setup is fresh /saline water interface at -50m, well from -45 to -42m. I want to see how the saline water will infiltrate the fresh water.

import numpy as np

import pygimli as pg
import pygimli.meshtools as mt
import pygimli.physics.ert as ert
import pygimli.physics.petro as petro
from pygimli.physics import ERTManager

################################################################################
# Create geometry definition for the modeling domain
world = mt.createWorld(start=[-20, 0], end=[20, -100], layers=[-20, -80],
                   worldMarker=False)

# make a filter at set set boundary condition name
line1=mt.createLine(start=[-20, -45], end=[-20, -42],boundaryMarker=12)
# Merge geometrical entities
geom = world + line1
pg.show(geom, boundaryMarker=True)

# Create a mesh from the geometry definition
mesh = mt.createMesh(geom, quality=32, area=0.2, smooth=[1, 12])
pg.show(mesh)

# Map regions to hydraulic conductivity in $m/s$
kMap = [[1, 1e-8], [2, 5e-3], [3, 1e-4]]

# Map conductivity value per region to each cell in the given mesh
K = pg.solver.parseMapToCellArray(kMap, mesh)

pg.show(mesh, data=K, label='Hydraulic conductivity $K$ in m$/$s', cMin=1e-5,
    cMax=1e-2, logScale=True, grid=True)

################################################################################
# The problem further boundary conditions of the hydraulic potential. We use
# :math:`p=p_0=0.75` m on the left and :math:`p=0` on the right boundary of the modelling
# domain, equaling a hydraulic gradient of 1.75%.

# Dirichlet conditions for hydraulic potential
left = 0.0
right = 0.75

pBound = {"Dirichlet": {12: left, 5: right, 6: right, 7: right}}

# Solve for hydraulic potential
p = pg.solver.solveFiniteElements(mesh, a=K, bc=pBound)

# Solve velocity as gradient of hydraulic potential
vel = -pg.solver.grad(mesh, p) * np.asarray([K, K, K]).T

ax, _ = pg.show(mesh, data=pg.abs(vel) * 1000, cMin=0.05, cMax=0.15,
            label='Velocity $v$ in mm$/$s', cMap='YlGnBu', hold=True)
ax, _ = pg.show(mesh, data=vel, ax=ax, color='k', linewidth=0.8, dropTol=1e-5,
            hold=True)

S = np.zeros((mesh.cellCount(),1))
# Fill injection source vector for a fixed injection position
ix=np.asarray(mesh.cellCenters())

# Below -50m, is the saline water
sourceCell=np.where(ix[:,1]<-50)[0]
S[sourceCell]=1  # I do not understand the unites here/
S=pg.Vector(S)

# Choose 800 time steps for 6 days in seconds
t = pg.utils.grange(0, 6 * 24 * 3600, n=800)

# Create dispersitivity, depending on the absolute velocity
dispersion = pg.abs(vel) * 1e-2

# Solve for injection time, but we need velocities on cell nodes
veln = mt.cellDataToNodeData(mesh, vel)
c = pg.solver.solveFiniteVolume(mesh, a=dispersion, f=S, vel=veln, times=t,
                             scheme='PS', verbose=0)

# We can now visualize the result:
i=0
for ci in c[1:][::20]:
    pg.show(mesh, data=ci * 0.001, cMin=0, cMax=3, cMap="magma_r",
        label="Concentration c in $g/l$",title=t,savefig='%i.png'%i)

    i=i+1

`

I started from your tutorial, but I am not sure if I am on the right path. My line of thought is that if Neumann BC should be applied above and bellow the filter. Any comments/suggestions?

[halbmy]

Dear Marios, very interesting issue. Probably boring for experienced groundwater modellers but not for us who know the same equation from DC modelling but are not so experienced in the details (BC and source) of the hydraulic and transport problems.

As for the hydraulic / Darcy problem, I think that's correct. Personally I would rather put a negative value on the filter screen and zero values on the right model boundaries, but the solution will be the same. All other boundaries (internal and external ones) are homogeneous Neumann by nature and that's ok in my opinion. Particularly, for the boundaries 1, 2, and 3 this makes sense as there is the symmetry axis and the water will flow towards the filter screen.

ax, _ = pg.show(mesh, data=pg.abs(vel) * 1000, cMin=1e-3, cMax=0.1, logScale=0,
                label='Velocity $v$ in mm$/$s', cMap='YlGnBu', hold=True,
                orientation='vertical')
pg.viewer.mpl.drawStreamLines(ax, mesh, p, color='k')

The streamlines look good to me.

[mariosgeo]

Thanks Thomas for your reply. Could you explain the units in the source term?

S[sourceCell]=1 

I did not divide by the area of the element, but is (1/area) actual 1g/l of something?

What if I want to have only one element that represents a heating wire? Then the

sourceCell = mesh.findCell([x, y])
S[sourceCell.id()]=1 

Corresponds to what units? Any hints how I could simulate heating for 10min this cell, and simulate an increase in temperature by dT?

I understand that pygimli is not a a primer hydrogeological tool, still it would be nice to have some basic modeling in this suite. In case, thanks a lot for your help and this beautiful package. It means a lot to our community.

[halbmy]

Now to the second part: You have a source vector that could be more easily created by

S = pg.Vector(mesh.cellCount())
S[pg.y(mesh.cellCenters()) < -50] = 1  # I do not understand the unites here/
pg.show(mesh, S)

Actually, you could also specify this as initial values for the concentration using u0=S instead (then the value would make sense as concentration unit), but I learned that it is also common to do this with the source term, i.e. additional saltwater flowing into the (initial) saltwater layer.

I've tested both variants (using S as f and u0). Both the u0 variant and the Svariant look ok to me, but are slightly different, not to speak about the unit tweaking of the latter.

Note that in general, the concentration affects density and therefore also the flow model, so one would have to iterate step by step over both hydraulic and transport equations. I guess that the influence will be limited in this scale, but of course in a larger scale, i.e. when modelling a freshwater lens.

[halbmy]

Well, according to the formula in https://www.pygimli.org/pygimliapi/_generated/pygimli.solver.html?highlight=finitevolume#pygimli.solver.solveFiniteVolume the source term f represents a change in u (here concentration) per time interval, so actually we should rather specify df/dt in order to be independent on the temporal discretization.

Note that in the S variant, we add more and more saltwater to the system which is responsible for the large numbers. So I guess the final state is visually ok, but not really (we should sum up the source terms and normalize). In contrast, in the u0 variant, the saltwater is slightly diluting as its area increases, which is also not fully correct.

Finally, what you're trying to model is actually a 3D problem that should be solved either in 3D, or in cylindersymmetric meshes, or by wavenumber transform (just the way we do it for the DC problem).

[halbmy]

Any more news on this issue? Otherwise we can close it (and keep is as a base for further discussions).

[halbmy]

Thank's for the interesting discussion, it would be interesting to compare the results to some "standard software". If this works out, let's try to generate an example for it.

Closing due to inactivity, but we can continue to discuss here.